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ABSTRACT 


A set of computer programs has been developed for the calcu- 
lation of the geomagnetic field and the tracing of field lines in 
space. The basic subroutine, geocentric ALLMAG, contains 
coefficients for seven recently-published field models as built- 
in data statements. At execution time the user can vary the 
model and/or the time period simply by changing input pa- 
rameters. Subroutine GDALMG is adapted for input and output 
in geodetic coordinates. ALLMAG and GDALMG are equivalent 
to Cain's FIELD and FIELDG, with the added flexibility of the 
choice of seven models. LINTRA traces field lines from any 
point in space to a specified altitude intersect in the same or 
opposite hemisphere, using any of the models contained in 
ALLMAG. Input is in either geocentric or geodetic coordi- 
nates, and output is returned in both. Mcllwain's INVAR 
package, which calculates B and L, has been adapted to use 
ALLMAG. All programs are described in detail, and sample 
calculations are given. 
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ALLMAG, GDALMG, LINTRA: COMPUTER PROGRAMS FOR 

GEOMAGNETIC FIELD AND FIELD-LINE CALCULATIONS 


1. INTRODUCTION 

The proliferation of quantitative geomagnetic field models in the last decade, 
evidently stimulated by the advent of the space age and the satellite era, has 
resulted in a confusing abundance of good numerical models that are based on 
a spherical harmonic expansion of the geomagnetic potential. Most of these 
models have included first- and occasionally second-order time derivatives of 
the spherical harmonic coefficients, giving the secular change of the field. 

There has been for some time a need for a versatile, unified set of computer 
routines which would permit the user to choose different models and/or time 
periods during one execution time, in order to compare the predictions of 
different models or to compute the value of the field at different time periods. 

In the past, when it was desired to perform field calculations with a different 
model or another time period, it was in many cases necessary (e.g. Mcllwain's 
INVAR) either to replace completely the field-computing routine or to first 
calculate the coefficients for the new date, punch them onto cards, and then 
physically insert them into the model deck; this was at best an error-prone and 
time-consuming procedure, lacking both flexibility and efficiency. 

The purpose of our effort was to devise a practical and dependable system by 
which these shortcomings could be overcome without sacrificing either accuracy 
or speed. It resulted in the creation of ALLMAG, a program package combining 
in one operation two previously separate functions, namely the calculation of 
coefficients for a given model and time period and the computation of the field 
vector. ALLMAG has the added advantage of incorporating seven widely used 
field models "under one roof," without having to read in data cards at execution 
time. These models, all with internal source terms only, are described in the 
next section. During a single execution one can successively vary either the 
model or time period, or both. 

The spherical harmonic coefficients for the models are stored in data statements 
in the beginning of the program; in this respect, the code is complete and self- 
contained. The coefficients are tested for accuracy and renormalized the first 
time the subroutine is called; updated coefficients are then calculated each time 
a new model or time period is selected. 
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ALLMAG is available in a "short version" (using subscripted variables and 
do-loops) and a time-saving "long version," with fixed indices and functions 
written out in extenso, yielding identical results. The long version calculates 
the field vector in about a third the time as the short version, and is therefore 
suitable for lengthy calculations with frequent calls to the subroutine. Execu- 
tion times are comparable to or faster than other programs currently available. 

ALLMAG was designed as an independent, self-contained subroutine with input 
and output in geocentric coordinates. A short subroutine GDALMG was de- 
veloped, with input and output in geodetic coordinates (geodetic latitude and 
altitude above the geoid) , to be used in conjunction with ALLMAG. GDALMG 
performs all necessary transformations between geodetic and geocentric coor- 
dinates and computes declination and inclination in addition to the geodetic field 
components. 

ALLMAG and GDALMG perform basically the same calculations as the MAGNET 
subroutine of Mcllwain's INVAR (McDwain 1961), as Hassitt and Mcllwain's 
NEWMAG (Has sitt and Mcllwain 1967 ), and as Cain's FIELD and FIELDG ( Cain 
et al . 1968), with the added advantage of offering a choice of seven models and 
different time periods at execution time. 

Two modified versions of ALLMAG were further developed, designed to be 
substituted for ALLMAG for special purposes. The first, DEKMAG, contains 
no built-in coefficients, but reads up to seven sets of coefficients as data cards 
at execution time. DEKMAG can thus accommodate new models as they become 
available. The second, ONEMAG, contains built-in coefficients for only one 
model, namely the IGRF model, and is designed for those who do not need the 
variety of models provided by ALLMAG. 

Finally, for line-of-force integrations and for computations of field-line inter- 
sects at any altitude level, a new field-line-tracing routine LINTRA was de- 
veloped, greatly improving Stassinopoulos 1 (1968) old LINTRA code. Designed 
around ALLMAG, the new LINTRA accepts initial starting positions in either 
geocentric or geodetic coordinates, and returns the initial position as well as 
the conjugate intersect in both systems. 

Special features of LINTRA include: (1) input control of tracing direction, 
permitting field-line tracing to the opposite hemisphere or to a lower altitude 
in the same hemisphere; (2) computation of arc length of line segment traversed; 
(3) search for and retention of minimum-B equator position and magnitude, if it 
has been crossed; (4) internal functional determination of optimum integration 
step size. 
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All programs are written in FORTRAN IV computer language and card decks 
are available in either the 029 model IBM keypunch format (EBCDIC) for use 
with the IBM 360 series machines, or the 026 keypunch format (BCD) for other 
FORTRAN-compatible computers. An INVAR program ( Hassi t t and M cI Iwain , 
1967) has been modified slightly so as to accept ALLMAG in place of NEWMAG. 

It yields identical values of B and L if used with the same model and time period. 
Sample calculations are given for GDALMG- ALLMAG, LINTRA, and the modi- 
fied INVAR. 

Complete listings are given for all programs except for INVAR (see Hassitt and 
McIIwain , 1967, for listings) and the long version of ALLMAG (see Cain et al . , 
1968, for listings of the arithmetic statements of the long version). 


2. SELECTION OF MODELS 

Seven geomagnetic field models were selected for inclusion into ALLMAG. 

Table 1 lists these models together with some of their basic characteristics. 

We included only models for which documentation was readily available through 
journal publication or otherwise; no preliminary models were included. We 
concentrated on the most recently-published models, although a few earlier, 
widely-distributed models were also included so as to enable direct comparisons 
with earlier subroutines. We did not include the widely-quoted Jensen and Cain 
model, since it contains no time- derivative terms and gives extremely poor 
predictions of the present-day field (see C ain and Sweeney , 1970). 

Note that Table 1 gives both the epoch and the data range. The term "data range" 
is used here to denote the time period during which geomagnetic data were ob- 
tained to define the model. The "epoch" of a model with secular time-derivative 
terms is the zero of time from which At is calculated in order to add or subtract 
the time-derivative contribution to the main-field terms. Since in this sense the 
epoch is simply a numerical constant, it may or may not lie within the data 
range. All the Cain models, for example, are based on an epoch of 1960, for 
simplicity, even though the recent POGO models utilize data only from POGO 
satellites, the first of which was launched in 1965. 

It is customary to set the input variable TM equal to the time period for which 
one wishes to calculate the field. It is generally undesirable, however, to input 
a time more than a few years away from the data range of a given model, since 
to do so requires extrapolation well outside the time period over which data 
was obtained to define the model. Recent studies ( Mead , 1972) have shown that 
such large extrapolations can lead to highly unreliable and often divergent re- 
sults. If we desire to predict the characteristics of the field in 1973.0, for 
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Table 1 

Geomagnetic Field Models La ALLMAG 


Model 

No. 

Designation 

Epoch 

Data Range 

n 1 

max 

No. of 
Coefficients 

— 

Reference 

Earth 

Documentation 2 

1 

GSFC 9/65 

1960 

1945-1964 

9 

99 

Oblate 

Hendricks & Cain, 1966 

2 

GSFC 12/66 

1960 

1900-1965.8 

10 

120 

Oblate 

Cain et al., 1967 

3 

POGO 10/68 

1960 

1965.8-1967.9 

ll 

143 

Oblate 

Cain and Langel, 1968 

4 

POGO 8/69 

1960 

i 

1965.8-1968.4 

10 

120 

Oblate 

Cain and Sweeney, 1970 

5 

IGRF 1965.0 

1965 


8 

i 

80 

Oblate 

IAGA Commission 2, 1969 
Cain and Cain, 1971 

6 

LME 1965 

1965 

1945-1964.5 

8 

80 

Spherical 

Leaton et al., 1965 

7 

USC&GS 1970 

1970 

1939-1968 

12 

168 

Oblate 

Hurwitz, 1970 
Hurwitz & Fabiano, 1969 


1 The FORTRAN-stored NMAX = n + 1, to avoid zero indices. 

max 7 

2 See References. 
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example, using a model such as POGO 10/68, whose data range is 1965.8-1967.9, 
it would be better to choose something like 1969.0 as a time input to the model, 
rather than 1973.0, which would require linear extrapolation over five years 
outside its data range. There are two reasons why these large extrapolations 
are undesirable. First, most models assume that the secular variation is linear, 
whereas long-term studies (Cain and Hendricks , 1968) clearly show that the 
secular variations are often highly non-linear and that linear extrapolations 
many years into the future are unreliable. Secondly, several of the models have 
used a relatively short time period to determine the secular time derivatives. 

It appears that data ranges of at least five years are necessary to clearly es- 
tablish the linear trends. 

Since the models contained in ALLMAG include no external sources, they yield 
very unreliable results beyond geocentric distances of 3-5 earth radii. Strong 
perturbations (10-100% or more) of the geomagnetic field are present in the 
outer magnetosphere. These perturbations depend on local time and season as 
well as solar wind conditions. Some improvement in the accuracy of field pre- 
dictions at these distances can be obtained by adding the contribution of external 
sources predicted by a model such as Mead's (1964), An improved model of the 
external field based on least-squares fits to satellite magnetometer data ( Mead 
and Fairfield , 1971), explicitly incorporating seasonal effects caused by the 
varying tilt of the dipole with respect to the solar wind, will soon be available. 

A few comments on each model are appropriate. Model 1 (GSFC 9/65) was based 
mostly on surface survey data, with some additional localized satellite data from 
Vanguard 3 and Alouette. This model, updated to the time period 1965, was in- 
corporated into the NEWMAG subroutine of the INVAR B-L program (Hassitt and 
McHwain, 1967). Model 2 (GSFC 12/66) incorporates survey data back to 1900 
and is the only model containing second time-derivative terms, thus making a 
quadratic fit to the secular variation. As such, it is probably the best single 
model fitting surface data over the period 1900-1965 (C ain and Hendricks , 1968). 
However, more recent models give better predictions of the current field. 

Model 3 (POGO 10/68) was the first model to use only satellite data (measure- 
ments of field magnitude only; no directional data included). Its data range is 
very short (2.1 years), and therefore its time derivative terms are rather 
poorly determined. With this model, the variable TM should probably be limited 
to the time period 1964-1969. Model 4 (POGO 8/69) is the most recently- 
published model from the Cain group as of this writing. Model 5 (IGRF) is now 
internationally-accepted as a reference field to use as a standard whenever 
comparisons are needed. Since it was derived from the components of many 
different models, no data range for it is given in Table 1. Model 6 (LME 1965) 
was used for the preparation of world magnetic charts for the epoch 1965.0 pub- 
lished by the Hydrographic Office of the British Ministry of Defence. The 
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British update their charts every 10 years, and a new model is not expected 
until 1975. The LME model is the only one included in ALLMAG whose deriva- 
tion neglected the oblateness of the earth (it assumed R = 6371.2 km every- 
where on the surface). Model 7 (USC&GS 1970) is the American World Chart 
Model for 1970, and was used to prepare the magnetic charts issued by the U.S. 
Naval Oceanographic Office in 1970. Models 6 and 7, based almost entirely on 
ground data, are probably preferable to use for predictions of the field at or 
near the Earth's surface. Models 3-5 would probably be preferable to use in 
space applications. Cain (1971) has recently summarized the status of geo- 
magnetic models obtained from satellite surveys. 


3. APPLICATION AND USE 

In this section the five basic routines are briefly discussed, their operation is 
described, and the arguments of their transfer vectors are presented. A more 
detailed description of the three essential programs ALLMAG, GDALMG, and 
LINTRA is given in the Appendix, with an analysis of the methods employed 
and a review of the organization and structure of the codes. Program listings 
and sample outputs are given in the attachment. 

A. ALLMAG 


Subroutine ALLMAG is the most essential part of the entire program package: 
it calculates the spherical geocentric field vector components and the total field 
strength at any position defined in geocentric coordinates, when given a time, a 
model number, the geocentric distance to the position, and the values of four 
trigonometric functions derived from its geocentric latitude and longitude. 

The variables in the calling sequence are 

ALLMAG(MODEL, TM, RKM, ST, CT, SPH, CPH, BR, BT, BP, B) 

The first seven arguments are input data: 

MODEL : integer from 1 to 7 which chooses the field model (see 
Section 2); 

TM : time for which the field coefficients are to be calculated, in 

decimal years A.D. (e.g., 1971.2); 

RKM : geocentric distance to given position, in kilometers; 
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The last four arguments are output data: 

BR, BT, BP : geocentric spherical field vector components (B r , positive 
outward; B 0 , positive southward; and B^, positive eastward) 
at the given position, in gauss; 



B : total magnetic field strength at the given position, in gauss. 

All arguments except MODEL are floating point variables. There are no re- 
strictions or limitations on the values of the input arguments other than on 8 , 
which at the poles (9 = 0,180) results in a divide by zero, and MODEL, which 
must be an integer from 1 to 7. 

For greatest efficiency, all calculations should be completed with one model and 
one time before changing either or both. 

There are no READ or other input statements in ALLMAG. A PRINT statement 
is executed each time a new model or new date is inputted. In addition, a PRINT 
statement is executed if any of the coefficients fail the error test in the initiali- 
zation process. 

Figure 1 shows a flow diagram of the program. A complete listing of the short 
version is given in the attachment. The cards of each deck are appropriately 
labeled in columns 73-80 as ALMGSxxx (short version) or ALMGLxxx (long ver- 
sion) where the three last columns (xxx) contain the sequential numbering, which 
for the short version is incremented by 2 but for the long version by 1, in order 
to accommodate the 635 cards. 

B. DEKMAG 

Subroutine DEKMAG is designed as a substitute for ALLMAG in order to accom- 
modate new geomagnetic field models as they become available. Instead of 
having coefficients built-in as DATA statements, it reads in coefficients for up 
to seven models as card input in standard format at execution time. Input decks 
corresponding to three of the ALLMAG models (IGRF, POGO 10/68, and POGO 
8/69) are included in the package. The calling sequence is identical to ALLMAG, 
with each variable having the same meaning. 
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When DEKMAG is called as a subroutine for the first time, a READ sequence 
is initiated. A deck containing one parameter card plus from one to seven sets 
of coefficients are to be placed in the card input stream so as to be available 
when the subroutine is first called. The input stream is as follows: 

Card 1: NMODLS (Col. 5), NPRINT (Col. 10) 

Deck A: First model to be read in 

Deck B: Second model, etc. 

The variable NMODLS, an integer from 1 to 7, controls the number of models 
to be read in. NPRINT controls the amount of information on each model printed 
in the output stream. If NPRINT = 0, only the label, as read in on the first card 
of the model deck, and the value of G(2,l), the largest coefficient, is printed out. 

If NPRINT = 1, the values of all coefficients are written on the output stream. 

The decks containing the model coefficients should be in the standard format as 
distributed by the Cain group at GSFC. The format for each model deck is as 
follows (see also Cain et al ., 1968): 

First card: Epoch (Cols. 4-9), Label (Cols. 10-73) 

Intermediate cards: N (Cols. 1-3), M (Cols. 4-6), GNM, HNM, GTNM, HTNM, 

GTTNM, HTTNM (11 Cols, each) 

Last card: Zero or blank in Cols. 1-3 

The values of the parameters J and K, normally punched in Cols. 1 and 2, re- 
spectively, of the first card of the model decks distributed by the Cain group, 
are assumed by the program to be zero; thus the program assumes that the co- 
efficients are for an oblate earth and are Schmidt-normalized. One should not, 
therefore, read in the Jensen and Cain coefficients (K = 1) or the Leaton, Malin, 
and Evans coefficients (J = 1) without modifying the program. Likewise, if 
DEKMAG is to be used with GDALMG, the branching statements in GDALMG 
for MODEL = 6 should be modified, unless the Leaton et al. model is to be read 
in as the sixth model. 

The models are assigned a number corresponding to the order in which they 
appear in the input stream. The integer MODEL in the calling sequence then 
determines which model is to be used. A STOP command is encountered if 
MODEL < 1 or MODEL > NMODLS. Each time a new model or new time is 
selected, a statement is printed in the output stream with the model number, 
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label (taken from the first card of each model deck), and time selected. For 
greatest efficiency, all calculations with one model and one time should be com- 
pleted before changing model or time. 

DEKMAG is supplied in the short version (with subscripted variables and do 
loops) and single precision. A long version can be assembled by replacing cards 
DEKMG 200-270 with cards ALMGL 222-635. Instructions are given in the 
listings for converting to double precision. 

C. O NE MAG 

Subroutine ONEMAG is designed as a substitute for ALLMAG where, for sim- 
plicity, only one model is required. The coefficients for the IGRF 1965.0 model 
are built into the routine as DATA statements. There are no READ statements 
in the program. A model other than the IGRF can be substituted by selecting 
the appropriate DATA statements from ALLMAG and relabeling the coefficient 
matrices as LG and LGT. 

The calling sequence is identical to ALLMAG except that the variable MODEL 
is removed. The routine is supplied as the short version, single precision. 

D. GDALMG 

Subroutine GDALMG represents a geodetic version of the magnetic field program. 
It actually is a short preparatory code that functions as an interface between a 
geodetic frame of reference and geocentric ALLMAG. Thus, it converts geo- 
detic latitude and altitude above the geoid into geocentric colatitude and geo- 
centric distance to the position and calculates the trigonometric functions for 
input into ALLMAG. The returning geocentric field vector components are 
transformed into geodetic components and the declination, inclination, and total 
horizontal field intensity are computed. The reference geoid is that adopted by 
the I.A.U. in 1964. A time and a model number have to be supplied to GDALMG 
for transmittal to ALLMAG. 

The variables in the calling sequence of this routine are 

GDALMG (MODEL, TM, GDLAT, GDLON, GDALT, X, Y, Z, F, H, DEC.AINC) 
where the first five arguments are input data: 

MODEL : same as in ALLMAG; 

TM : same as in ALLMAG; 
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Figure 2. Flow Diagram of GDALMG 
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GDLAT : geodetic latitude, in degrees; 

GDLON : east longitude, in degrees (invariant in the two coordinate 

systems); 

GDALT : geodetic altitude or height above geoid, in kilometers. 

The last seven arguments are output data: 


•X, Y, Z 

F 

H 

DEC 

AINC 


geodetic north (%-B#), east (=130), and downward vertical (%-B r ) 
field components at the given position, in gauss; 

scalar magnitude or intensity of field at the given point, in gauss; 

total horizontal field intensity, in gauss (H 2 = X 2 + Y 2 ); 

declination, in degrees (-180 < DEC < 180), positive eastward; 

inclination or dip angle, in degrees (-90 < AINC < 90), positive 
downward. 


As in ALLMAG, all arguments except MODEL are floating point variables which 
have no restrictions on the values they may assume, except that GDLAT = ±90° 
will result in a divide by zero, and MODEL must be an integer from 1 to 7. As 
with ALLMAG, for greatest efficiency all calculations should be completed with 
one model and one time before changing either or both. There are no input- 
output statements in GDALMG. See ALLMAG section, however, for a descrip- 
tion of its output statements. 


A flow diagram of GDALMG is given in Figure 2 and a listing of the code is 
included in the attachment. The deck is labeled in Columns 73-80 as GDALMxxx, 
with the last three columns containing a sequential numbering of the cards by 
increments of 2. 


E. LINTRA 

The LINTRA code is an independent "Main" program that traces a field line 
passing through a given point on or above the geoid along either direction. In 
this process the program obtains intersects at any specified altitude level(s), it 
calculates and retains field strength and position of the minimum-B equator (if 
crossed), and it computes the arclength of the line- segment traversed. Designed 
around ALLMAG, the code accepts input coordinates for starting positions in 
either a geocentric or a geodetic system and it returns the coordinates for the 
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Figure 3. Flow Diagram of LINTRA 
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desired intersect in both systems. The tracing direction is input controlled; 
for a space point, this permits tracing the field line to the opposite hemisphere 
or down the same hemisphere, or both. The code requires a constant integra- 
tion step size during execution time, but the increment may be changed from 
run to run. For maximum flexibility and efficiency the user has the option of 
an optimum step size calculated internally as a function of magnetic latitude 
and geocentric distance or of determining his own step size through input. In 
this case, lines 132 and 134 of the program should be commented out because 
in its present form, the code is designed to override the given input step size. 

The code offers two output choices: either a continuous step by step account of 
the tracing process, where the coordinates of each successive point are printed 
• with the positional total field strength and its components, plus the final results, 

or the final results only. 

LINTRA has two supporting subroutines: 

1. ITERAT, which performs the actual integration of the differential 
equations of the field line after a 7- step initialization process during 
which the tracing point is advanced by three increments, and 

2. CONVRT, which performs all programmed conversions from geocentric 
to geodetic coordinates or vice versa. 

The user-supplied input data for LINTRA are: 

First input card (FORMAT (15, F10.2, 215) ): 

MODEL : same as in ALLMAG; 

TM : same as in ALLMAG; 

NPRINT : output control: 

= 1 : prints each integration step, plus final results; 

^ 1 : no r unnin g print-out, only final results; 

ICOORD : reference system of input coordinates: 

= 1 : geodetic; 

= 2 : geocentric; 

Each successive card (FORMAT(6F10.6, 2A4) ): 
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If ICOORD = 1: 


GDLAT, GLON, ALT 


geodetic latitude and east longitude, in degrees, 
and altitude above geoid, in kilometers; 


Or, if ICOORD = 2: 


GCLAT, GLON, RKM : 


geocentric latitude and east longitude, in degrees, 
and geocentric distance, in kilometers; 


DS 


integration step size (tracing increment), in 
kilometers; 


(Note: This parameter may be omitted, since the program determines DS unless 
cards LINTR 132 and 134 are commented out.) 


tracing-direction control: 

= >0. : traces towards higher altitudes; 

= <0. : traces towards lower altitudes; 

geodetic altitude of desired intersect, in kilometers; 

name of station or designation of origin (starting 
point). 

MODEL, NPRINT, and ICOORD are integers; LABEL1 AND LABEL2 are alpha- 
numeric; all other arguments are floating-point variables. If DIR = -Land 
HALT > ALT, the line is not traced. 


DIR 

HALT 

LABEL1, LABEL2 


If long output has been specified (NPRINT = 1), the code will print the following 
variables for every integration step: 


L 

DLATP, DLONP 

RP, HP 

BR, BT, BP 
B 


: current step- count; 

: geocentric latitude and longitude of last position 

(L-l), before iteration, in degrees; 

: geocentric distance and altitude of last position 

(L-l), before iteration, in kilometers; 

: same as in ALLMAG, at last position (L-l); 

: same as in ALLMAG, at last position (L-l). 
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Line-tracing proceeds until either (a) the number of steps exceeds 200 or (b) the 
intersect is reached on the downward path (i.e., for R < RP because the tracing 
skips over the intersect in the same hemisphere, where R > RP). If a HALT > 
initial ALT in the same hemisphere is desired, card LINTR200 should be re- 
moved or commented out. 

When the specified intersect has been reached, the final results are printed: 

PLAT, PLON, PRKM : geocentric coordinates of intersect, in degrees 

and kilometers; 

PGDLAT, PLON, PGALT : geodetic coordinates of intersect, in degrees 

and kilometers above the geoid; 

ARC : length of traced field line segment, in 

kilometers. 

A flow diagram of LINTRA is given in Figure 3 and a listing of the program is 
contained in the attachments. 


ACKNOWLEDGMENTS 

We are indebted to Joseph Cain, Robert Langel, and Ronald Sweeney for discus- 
sions of the characteristics of various models and for the use of their routines 
FIELD and FIELDG, in both the long and short versions, which served as the 
basic framework for ALLMAG. We thank Eugene Fabiano and Louis Hurwitz 
of the Coast and Geodetic Survey for their deck to compute the USC&GS Field 
and for their testing of ALLMAG on a CDC 6600. 

We further thank A1 Geelhaar for his very valuable programming assistance, 
and Carl Mcllwain and Gerald Peters of the University of California, San Diego, 
for testing ALLMAG on their computers. 


16 


REFERENCES 


Cain, Joseph C., Shirley J. Hendricks, Robert A. Langel, and William V. 

Hudson, A proposed model for the International Geomagnetic Reference 
Field, 1965, J. Geoma gn. Geoelec ., 19, 335-355, 1967. 

Cain, Joseph C., and Shirley J. Hendricks, The geomagnetic secular variation 
1900-1965, NASA Technical Note TN D-4527, April 1968. 

Cain, Joseph C., Shirley Hendricks, Walter E. Daniels, and Duane C. Jensen, 
Computation of the main geomagnetic field from spherical harmonic ex- 
pansions, Data Users' Note NSSDC 68-11, National Space Science Data 
Center, Greenbelt, Maryland, May 1968. 

Cain, Joseph C., and Robert A. Langel, The geomagnetic survey by the Polar 
Orbiting Geophysical Observatories OGO-2 and OGO-4 1965-1967, GSFC 
Report X-612-68-502, Greenbelt, Maryland, 1968. 

Cain, Joseph C., and Ronald E. Sweeney, Magnetic field mapping of the inner 
magnetosphere, J. Geophys. Res. , 75^ 4360-4362, 1970. 

Cain, Joseph C., Geomagnetic models from satellite surveys, Rev. Geophys . 
and Space Phys ., 9, 259-273, 1971. 

Cain, Joseph C., and Shirley J. Cain, Derivation of the International Geomag- 
netic Reference Field (IGRF 10/68), NASA Technical Note TN D-6237, 
August 1971. 

Chapman, S., and J. Bartels, Geomagnetism, Oxford University Press, London, 
1940. 

Hassitt, A., and C. E. Mcllwain, Computer programs for the computation of B 
and L (May 1966), Data Users' Note NSSDC 67-27, National Space Science 
Data Center, Greenbelt, Maryland, May 1967. 

Hendricks, S. J., and J. C. Cain, Magnetic field data for trapped-particle 
evaluations, J. Geophys. Res ., 71, 346-347, 1966. 

Hurwitz, Louis, Mathematical model of the 1970 geomagnetic field, ESSA Coast 
and Geodetic Survey, preprint, May 4, 1970. 

Hurwitz, L., and E. B. Fabiano, Geomagnetic secular variation 1937.5-1967.5, 
ESSA Coast and Geodetic Survey, Preprint, August, 1969. 


17 


IAGA Commission 2 Working Group 4 Analysis of the Geomagnetic Field, 

International Geomagnetic Reference Field 1965.0, J. Geophys. Res., 74, 
4407-4408, 1969. 

Leaton, B. R., S. R. C. Malin, and Margaret J. Evans, An analytical representa- 
tion of the estimated geomagnetic field and its secular change for the epoch 
1965.0, J. Geomagn. Geoelec. , 17 , 187-194, 1965. 

Mcllwain, C. E., Coordinates for mapping the distribution of magnetically 
trapped particles, J. Geophys. Res., 66, 3681, 1961. 

Mead, Gilbert D. , Deformation of the geomagnetic field by the solar wind, 

J. Geophys. Res., 69, 1181-1195, 1964. 

Mead, Gilbert D., and Donald H. Fairfield, Quantitative magnetospheric models 
derived from satellite magnetometer data, Trans. A.G.U., 52, 318, 1971. 

Mead, Gilbert D., Secular change of geomagnetic conjugate points, preprint, 
1972. 

Morrison, John, and Samuel Pines, The reduction from geocentric to geodetic 
coordinates, Astron. J., 66, 15-16, 1961. 

Roederer, J. G., W. N. Hess, and E. G. Stassinopoulos, Conjugate intersects to 
selected geophysical stations, GSFC Report X-642-65-182, April, 1965; 
also in Goddard Space Flight Center contributions to the COSPAR meeting - 
May 1965, NASA Technical Note TN D-3091, July 1966. 

Stassinopoulos, E. G., Computer codes for geomagnetic field line tracing and 
conjugate intersect calculations, GSFC Report X-642-68-429, November, 
1968. 


18 



APPENDIX A 


COMMENTS ON ALLMAG 


ALLMAG is the fundamental routine which calculates the three components of 
the vector field from a spherical harmonic expansion. The relevant mathe- 
matical formulation is readily available in, for example, Chapman and Bartels 
(1940) and Cain et al . (1968), and we shall not repeat it here. We have adhered 
to Cain's treatment in most respects. ALLMAG is essentially similar in its 
basic mathematics to Cain's subroutine FIELD, with the following differences 
in the organization of the program: 

1. Any one of seven models may be selected; successive selection of 
different models is possible. 

2. All coefficients for the seven models are built into ALLMAG as DATA 
statements. ALLMAG contains no READ statements. 

3. An internal sum-test is performed on all the coefficients to check for 
accuracy the first time the subroutine is called, and all coefficients 
are then converted from Schmidt to Gauss normalization. 

4. Different time periods may be selected if desired. 

5. A statement indicating the model and time period selected is written on 
the output stream the first time ALLMAG is called, and each succes- 
sive time that either MODEL or TM is changed. 

6. ALLMAG is completely self-contained if input and output in geocentric 
coordinates is desired. No other supporting routines are needed. 
(GDALMG should be used with ALLMAG if geodetic input and output is 
desired.) 

All coefficients are stored in the DATA statements in their Schmidt-normalized 
form, since this is the normalization used in all recently-published models 
(i.e., in Cain's notation, K = 0 for all models). All models except Model 6 
( Leaton et al ., 1965) were derived using an oblate earth (J = 0 in Cain's notation). 
ALLMAG assumes that the appropriate conversion of input quantities into geo- 
centric coordinates has already been performed. See the appendix section on 
GDALMG for a discussion of these conversions. 
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The variable TM used as input to ALLMAG may specify a time that is earlier 
or later than the epoch assigned to each model; the time derivatives or secular- 
variation terms are then used to update or predate the coefficients. However, 
it is best not to input a time more than a few years away from the data range of 
a given model, i.e., the actual time period during which geomagnetic data were 
used to define the model. See Section 2 on models for further discussion of 
this point. 

ALLMAG is designed to handle any sets of coefficients with up to 168 terms, 
that is, n max = 12 (Fortran NMAX = 13), but the "short version" can be easily 
modified to accept models with a larger n max . 

The coefficients in the data statements have been checked very carefully against 
the published coefficients for each model, and we believe that there are no 
errors. To further reduce the possibility of errors due to misplaced data cards, 
that is, to insure that the data cards with the coefficients are in proper sequence, 
a special test was included in the code, which during initialization compares the 
sum of the products of the model-coefficients times their sequence number with 
precalculated values placed into a data statement under the variable name 
LSUM. This test is performed only once for all models during the first call, in 
the initialization step. If a discrepancy is detected, a message will be printed, 
indicating which model is being questioned and giving the precalculated as well 
as the computed sums; the run is then terminated. 

Comment statements throughout ALLMAG highlight significant parts or functions 
of the routine while all input/output variables, constants, and parameters are 
defined and identified in comment statements at the beginning of each routine. 

ALLMAG was tested and has successfully run on IBM-360 computers (models 
40, 75, and 91) in the 029 key-punch version (EBCDIC). It was also tested with 
equal success on an IBM-7094, a CDC-6600, and a UNIVAC-1108 computer in 
the 026 key-punch version (BCD). To run the program on the GSFC IBM-7094 
we had to trim the code by three models to decrease its size because of limita- 
tions in the system's capabilities with regards to table- space. 

Timing of the long version of ALLMAG on the IBM 360/91 produced the following 
results: 

Model Average Time Per Call 

340 microsecs 


IGRF (n = 8) 

'max 1 

GSFC 12/66 (n max = 10) 


20 


450 microsecs 


These times axe averages, obtained from five separate rims of one thousand B 
calculations each. They include the time used for initialization, coefficient 
recalculation, and I/O interrupt (but not actual I/O time). The short version 
times are about a millisecond per call. 

The model coefficients are stored as integers rather than real variables. This 
has two advantages: we found that compilation time is significantly less, and 
the sum-test, using integer arithmetic, can be exact. The integers are con- 
verted to real variables and divided by the appropriate power of 10 (stored as 
G1 (1,1), etc.) at the same time as they are converted from Schmidt to Gauss 
normalization. Equivalence statements are used to conserve storage space. 

The EBCDIC short version is supplied in single precision and the EBCDIC long 
version is supplied in double precision. Only two card changes are necessary 
to convert either version from single to double precision or vice versa; these 
changes are indicated in the comment cards. The BCD decks are supplied in 
single precision only. 
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APPENDIX B 
COMMENTS ON GDALMG 


The subroutine GDALMG accepts input in geodetic coordinates, i.e., geodetic 
latitude and altitude above the geoid. It converts these quantities into geocentric 
co-latitude and geocentric distance for input into ALLMAG. (Geodetic and geo- 
centric longitude are equivalent.) The geocentric field components obtained 
from ALLMAG are then converted into geodetic components, and the horizontal 
component, declination, and inclination are calculated. 

♦ The geometry is shown in Figure 4, where the earth is depicted as a grossly 

exaggerated ellipsoid. The altitude h is measured along a plumb line perpen- 
dicular to the geoid; and the geodetic latitude X. is measured at the intersection 
of the plumb line with the equatorial plane. 

The reference geoid is that adopted by the International Astronomical Union in 
1964. The two parameters defining this geoid are 

a = 6378. 16 km 

a - b 1 

f a 298.25 

where a is the equatorial radius, b is the polar radius, and f _ is the flattening 
factor. From these parameters we can derive the following quantities: 

b = 6356.7747km 

a 2 /b 2 = 1/(1 - f) 2 = 1.00673966 

e 2 = a 2 /b 2 - 1 = 0.00673966 

where e is the so-called second eccentricity. The geocentric latitude at the 
surface of the geoid §_ is related to the geodetic latitude by 


tan /3 


b2 

tan A. 

a 2 


from which 
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sin /3 


= sin A. 


A 


sin 2 \ + — 
b 4 


cos 





The geocentric distance to the geoid, p, is given by 

p - a/(l + e 2 sin 2 /3)^ 


The coordinates (x,y) of the space point P are then given by 

x = p cos /3 + h cos k 
y “ p sin j3 + h sin k 

and thus the geocentric distance R and colatitude 9. are 

R = (x 2 + y 2 ) 54 

6 - cos -1 (y/R) 

The small angle 8 at the space point between the downward vertical and the 
geocentric direction is given by 

8 = e + \ - 90° 

which is positive for positive latitudes and negative for negative latitudes. Thus 

sin 8 = sin 6 sin A. - cos 6 cos ^ 


cos 8 - cos 9 sin \ + sin 9 cos k 


All of the above equations are exact for any altitude and an ellipsoid of any 
eccentricity. 

The geodetic field components X (northward), Y (eastward), and Z (downward 
vertical) are related to the geocentric components B , B e , and B^ approximately 
by 
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x 




-B e 
Y = B 0 
Z is - B r 

but to get the exact relationship we must rotate B r and B 0 through the small 
angle 5: 

X = - Bg cos S - B r sin S 
Z = B q sin S - B r cos S 

Finally, the horizontal component H, declination D, and inclination I_ are defined 
by 

H = (X 2 + Y *f 

sin D = Y/H 

cos D = X/H 

tan I = Z/H 

Model number 6 ( Leaton et al . , 1965) was derived from surface measurements 
without correcting for the obiateness of the earth (i.e., in Cain's notation, J = 1 
for this model). Therefore, to use this model properly, altitudes should be 
referenced to a sphere of radius 6371.2 km instead of the geoid. For this model 
GDALMG sets the geocentric latitude equal to the geodetic latitude and sets 
R = 6371.2 + h. 

Although our equations are somewhat simpler, the geodetic transformations 
outlined here are exactly equivalent to those given by Cain et al . (1968) in de- 
scribing the corresponding portions of FIELDG. (Note, however, that the ex- 
pression for FAC at the bottom of their page 5 should be squared; it is correct 
in FIELDG.) Identical results are obtained from GDALMG and FIELDG if the 
same model and time period are used. 


26 



APPENDIX C 


COMMENTS ON LIN TEA 


Field-line tracing routines are useful in a number of applications: for example, 
in locating the position of conjugate points, in determining mirror point locations 
for particles at a given position with a given pitch angle, in calculating the value 
of the second adiabatic invariant prior to calculating the Mcllwain L-parameter, 
in ascertaining the topology of the magnetosphere, etc. Conceptually, the 
process is straightforward: given a point in space and a magnetic field model, 
we wish to trace a line passing through that point which is everywhere parallel 
to the local field direction as determined by the model. The brute-force tech- 
nique is to approximate the line by a series of short, straight line segments; 
each segment is parallel to the local field direction only at one end. Then at 
the other end of the segment, another straight line segment is constructed 
parallel to the field at that point. The true curve is approximated by making 
the segments as short as is required to achieve the desired accuracy. Such a 
technique is usually time-consuming because of the necessity of keeping each 
segment short in order to approximate the true curve. 

A much more satisfactory approach is to solve the differential equations defining 
a field line. In cartesian coordinates these are 

dx/ds = B x /B 

dy/ds = B y /B 

dz/ds = B z /B 

where ds is the infini tesimal arc length along the line. 

In spherical coordinates these equations become 

dr/ds = B r /B 

dd/ds = (l/r) (B*/B) 

d<£/ds = (l/r sin 6) (B^/B) 

The LINTRA line-tracing routine integrates these equations with the aid of 
subroutine ITERAT, us ing an A dams 4-point integration formula. Since each 
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integration step takes into account the local curvature of the line, a much larger 
step size can be used without sacrificing accuracy. 

An early version of this routine was used by one of us to calculate conjugate- 
point positions (Roed erer et al ., 3.965). This program was later documented in 
a modified form as LINTRA (St assinopoulos , 1968). It should be noted that our 
technique differs in principle from that used in the LINES subroutine of Mcllwain's 
INVAR program to calculate B-L coordinates. LINES varies the step size along 
the field line as its curvature changes; LINTRA, on the other hand, uses a con- 
stant, predetermined step size. 

The present LINTRA is essentially an improved and faster version of the old 
code, adapted to use the new ALLMAG and to which some desirable features 
were added, as for example the coordinate conversion facility, the output choice, 
the internal functional determination of optimum integration step size and 
others. Retained in the program was the search for and the storage of the 
minimum-B location encountered in the integration process; the respective 
geocentric coordinates are easily accessible for output in cards LINTR 174- 
178. It should be noted, however, that the saved minimum-B value may not 
correspond to the true minimum-B equator. In fact, the minimum-B equator 
may not lie at all between the origin and the conjugate intersect. In that case 
the stored minimum-B location will lie very close to either one of these two 
end-points of the field line. If crossed in the tracing process, the location of 
the minimum-B equator is approximated along the field line with a maximum 
error of plus or minus one arc increment. 

An addition to the new LINTRA is the calculation of an optimum step size D S 
as an empirical function of the geomagnetic dipole coordinates of the starting 
point: 


DS = .06 R/cos 2 \ d - 370 km 


where R is the geocentric distance in kilometers and k^ is the dipole latitude, 
defined by 

sin \ d = cos 11.5° sin A. g + sin 11.5° cos k^ cos (0 g + 69°) 

where \ g and $6 g are the geocentric latitude and east longitude, and where the 
geocentric position of the north dipole is assumed to be at 11.5° colatitude, 

69° west longitude. Since the empirical formula gives very large values near 
the magnetic poles, the step size is limited to a maximum of 3000 kms. We 
have found that the step size given by this formula is large enough to trace to 
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the opposite hemisphere in less than 200 steps, yet small enough to insure high 
accuracy. We have found that if high accuracy is not required, D S can be multi- 
plied by factors of 2 to 5, thus saving significant amounts of computer time. If 
desired, the user can input the step size, removing or commenting out cards 
LINTR 130-134. 

If one wishes to trace a line initially towards the earth's surface (DIR < 0), a 
shorter step size is usually desirable. In this case the empirical formula is 
modified so that DS is no greater than (ALT-HALT)/20, where ALT is the 
initial altitude and HALT is the final desired altitude. 

The old concept of a user-controlled tracing direction has not been changed. 

However, the range of applicability of new LINTRA was expanded so as to include 
intersects in the opposite hemisphere which were formerly excluded, namely 
intersects that have a higher altitude than that of the origin. 

In its present form, the program bypasses an intersect in the same hemisphere 
that lies above the initial point (HALT > ALT). But this bypass can be removed 
by commenting out card LINTR 200 from the deck. If the bypass is removed, 
and if the desired intersect lies close to the initial point, especially when large 
L-values are involved (high latitudes), it is necessary to override the internally 
calculated step size and input a step size small enough for effective tracing, 
that is, no larger than 0.1 x (HALT- ALT). 

Tracing is terminated whenever the altitude of the point on the field line becomes 
equal to or less than HALT (the specified altitude for the intersect). When this 
point is reached, the latitude and longitude at the conjugate intersect altitude h x 
(= HALT) are determined by a quadratic interpolation formula. If \ , X_ 2 , and \ 
are successive values of the latitude at altitudes h x , h 2 , and h 3 , which bracket 
the desired intersect altitude h x , then the latitude at h. is given by 

f • \ . r M 

_ (h 2 -h 3 )(h i -h 2 )(h i - h 3 )\ 1 H-(h 3 -h 1 )(h i -h3)(h i -h 1 )\ 2 +(h 1 -h 2 )(h i -h 1 )(h.-h 2 )A.3 
^ ~ .. . (hj-hj) (h^ha) (h 2 -h 3 ) ... 

The interpolated longitude is given by a similar formula. 

Finally, in regard to use and results it should be noted that: 

1. Tracing accuracy improves only slightly with decreasing integration 
step size while running-time increases disproportionally; therefore, 
the calculated step sizes are a good compromise,. ... 
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The location of conjugate points or intersects will vary with every field 
model; the dispersion with different models may be substantial; a "true" 
position, in an absolute sense, cannot be established. More accurate 
results can be obtained: 

a. For low L-values (L < 3), by using that model whose components 
differ the least from actual surface measurements, over adjacent 
areas of the globe; 

b. For high L-values (3 < L < 6), by the inclusion of external source 
terms into the model representation so as to account for contribu- 
tions resulting from the ring current in the magnetosphere and the 
currents on the magnetospheric boundary, which contributions be- 
come significant at L-values greater than 4; 

c. By considering and accounting for: 

(1) diurnal effects (local time effects) 

(2) severe magnetic storm effects 

(3) seasonal effects 

(4) solar cycle effects 

3. The indiscriminate application of the time derivatives of a field model 
, to predict the expected secular variations too far in the future may 
introduce a larger error into the calculations than expected (see Sec- 
tion 2). 

Subroutine ITERAT is a slightly modified version of the old DPNV program 
( Stassinopoulos , 1968). It uses a bootstrap method of getting started along the 
field line; three steps of length DS are taken in the first seven iterations. The 
subroutine stores the derivatives of the spatial quantities from the preceding 
steps, and henceforth one step along the curved field line is taken each time 
ITERAT is called. _ 

The subroutine CONVRT is used to transform coordinates from geocentric to 
geodetic or vice versa. The geodetic to geocentric conversion is identical to 
the one used in GDALMG. The geocentric to geodetic conversion uses the for- 
mulas given by Morr iso n and Pines (1961). The two parts are compatible, i.e., 
a geodetic to geocentric to geodetic conversion returns the same initial values 
with six-figure accuracy. 
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APPENDIX D 


i- — - - 


USE OF ALLMAG WITH INVAR 


For the benefit of the numerous users of Mcllwain's INVAR program, which 
calculates the magnetic parameter L, we have replaced subroutine NEWMAG 
(corresponding to the old MAGNET) with ALLMAG and named the modified deck 
"INVARA" (for INVAR - ALLMAG). Although very minor changes were made in 
three subroutines only, the letter A was added to the end of the name of all pro- 
grams in the deck to differentiate it from the standard versions. The routines 
INVAR, START, and LINES were modified. Their input-output arguments were 
implemented to carry the input parameters "MODEL" and "TM" needed by 
ALLMAG, while statements number 28-36 and 64-74 were added to START and 
statements 112-122 to LINES, in order to produce the appropriate calling se- 
quence for ALLMAG. Finally, a main program was constructed to read sample 
input points in either geocentric or geodetic coordinates, convert to the opposite 
coordinate system, call INVARA, and print results. The subroutine CONVRT 
(see LINTRA appendix) was added to the package to perform the conversions 
required by the main program. 
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SUBROUTINE ALLMAG< MODEL ,TM , RKM , ST ,CT , SPH ,CPH , BR , BT , BP , B ) ALMGS002 

*** GEOCENTRIC VERSION OF GEOMAGNETIC FIELD ROUTINE ALMGS004 

*** SINGLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC, 029 PUNCH) ALMGS006 

*** SHORT DECK, USES SUBSCRIPTED VARIABLES AND DO LOOPS ALMGS008 

*** EXECUTION TIME PER CALL 3 TIMES GREATER THAN LONG DECK ALMGS010 

*** PROGRAM DESIGNED AND TESTED BY' E G ST ASS I NOPOULOS AND G D MEAD, ALMGS012 

*** CODE 641, NASA GODDARD SPACE FLT CTR, GREENBELT, MD 20771 ALMGS014 

***** INPUT: MODEL CHOICE OF 7 MODELS - SEE BELOW ALMGSO 1 6 

***** RKM GEOCENTRIC DISTANCE IN KILOMETERS ALMGS018 

***** TM TIME IN YEARS FOR DESIRED FIELD ALMGS020 

***** ST »CT SIN 6 COS OF GEOCENTRIC COLATITUDE ALMGS022 

***** SPH,CPH SIN £ COS OF EAST LONGITUDE ALMGS024 

***** OUTPUT: BR,BT,BP GEOCENTRIC FIELD COMPONENTS IN GAUSS ALMGS026 

***** B FIELD MAGNITUDE IN GAUSS ALMGS028 

***** NOTE: FOR GREATEST EFFICIENCY, COMPLETE ALL CALCULATIONS WITH ALMGS030 
ONE MODEL AND ONE TIME BEFORE CHANGING MODELS OR TIME. ALMGS032 
***** for DOUBLE PRECISION ADD THE FOLLOWING CARD ALMGS034 

IMPLICIT REAL*8< A-H,0— Z ) ALMGS036 

***** see END OF DECK FOR ONE MORE CHANGE ALMGS038 

REAL*8 LABEL (4,7) / • HENDR ICKSCC A I N 99-TERM GSFC 9/65 CAIN ET. AL. ALMGS040 

A 120-TERM GSFC 12/66 CAIN6LANGEL 143-TERM POGO 10/68 CAIN&SWEENEY ALMGS042 
B120-TERM POGO 8/69 I GRF 1965.0 80-TERM 10/68 LEATON MALIN EVALMGS044 


CANS 80-TERM 1965 HURWITZ US C£GS 168-TERM 1970*/ ALMGS046 

DIMENSION TO ( 7 ) , NMX( 7),ISUM(7,3),G(13,13) ALMGS048 

DATA T0/4*1960. ,2*1 965., 1970. /,NMX/ 10, 11,12,11,9,9,13/ ALMGS050 

INTEGER LSUMI 7,3) /- 1646 106,-1 795 169, -1865298, -1777057, -158472, ALMGS052 

A-156856, -2 191 704,-62661,-96778, -18 15 19, -835 55, -9 569, -9599, ALMGS054 

B-8593, 1,-10618, 5*1/ ALMGS056 

INTEGER*4 G1(13,13),GT1(13,13),GTT1(13,13),G2(13,13),GT2(13,13), ALMGSO 58 

1 GTT2( 13, 13) ,G3( 13,13),GT3< 13, 13 ) ,GTT3 ( 13, 13 ) , G4 ( 1 3, 13 ) , ALMGS060 

2 GT4( 13,13) ,GTT4( 13,13 ) ,G5( 13, 13) ,GT5( 13, 13) ,GTT5( 13, 13) , ALMGS062 

3 G6( 13,13), GT6( 13,13) ,GTT6< 13, 13) »G7 ( 13,13 ) »GT7 ( 13, 13) »GTT7( 13, 13JALMGS064 

4 »LG< 13,13,7) ,LGT( 13,13,7) ,LGTT( 13, 13,7 ) ALMGS066 

REAL*4 GG< 13, 13,7), GGT ( 13,13,7 ),GGTT( 13, 13,7) ,SHMIT( 13,13) ALMGS068 

EQUIVALENCE ( G1 ( 1 > ,GG( 1 ) ,LG( 1 ) ) , ( GT1 ( 1 ) ,GGT ( 1 ) , LGT ( 1 ) ) , ALMGS070 


(GTTK 1 ) , GGTT ( 1 ) ,LGTT ( 1 ) ) , 


<G2< 1 ) »LG( 1,1,2)), 
(G3<1), LG< 1,1,3) ), 
( G4( 1 ) ,LG< 1,1,4)), 
( G5 ( 1 ) » LG( 1,1,5) ), 
<G6( 1 ) ,LG( 1,1,6)), 
( G7( 1 ) » LG( 1 » 1 » 7 ) ), 


<GT2< 1 ) , LGT ( 1,1,2)), 
( GT3 ( 1 ) , LGT ( 1,1,3)), 
(GT4( 1) ,LGT( 1,1,4) ) , 
<GT5(1)*LGT (1,1,5) ), 
( GT6 ( 1 ) , LGT ( 1,1,6) ) , 
( GT7 (1), LGT (1,1, 7) ), 


***** the FOLLOWING DATA CARDS CONTAIN THE FI 
***** FOR THE FOLLOWING SEVEN MODELS: 

***** G1 »GT1 : HENDRICKS £ CAIN 99-TERM 

***** G2,GT2,GTT2: CAIN ET. AL. 120-TERM 

***** G3»GT3 : CAIN 6 LANGEL 143-TERM 

***** G4,GT4: CAIN £ SWEENEY 120-TERM 

***** G5 »GT5 : I GRF 1965.0 80-TERM 

***** G6,GT6: LEATON MALIN £ EVANS 1965 

***** FOR MODEL 6 ( LME 1965) SET RKM 

***** G7,GT7: HURWITZ US COAST £ GEODETIC S. 


( GTT2 ( 1 ) , LGTT ( 1,1,2) ) , 
(GTT3(1 >,LGTT( 1, 1,3) ) , 
( GTT4 ( 1 ) , LGTT ( 1,1,4) ) , 
( GTT5 ( 1 ) , LGTT (1,1,5) > , 
( GTT6 ( 1 ) , LGTT I 1,1,6) ) , 
(GTT7 ( 1 ) , LGTT (1,1,7)) 
FIELD COEFFICIENTS 


***** G6,GT6: LEATON MALIN £ EVANS 1965 80-TERM EPOCH 196 

***** FOR MODEL 6 (LME 1965) SET RKM = 6371.2 + ALTITUDE 

***** G7,GT7: HURWITZ US COAST £ GEODETIC S. 168-TERM EPOCH 197 

DATA G1 / 10, -304249,-15361,13009,9576,-2277,498,709,48,99,3*0, 
A 57748,-21616, 30002 ,—19870, 8028,3595,6 07, -572, 67, 29, 3*0, -19498, 


GSFC 9/65 
GSFC 12/66 
POGO 10/68 
POGO 8/69 
10/68 
80-TERM 
= 6371.2 + 
168-TERM 


T ( 13,13) ALMGS068 

), ALMGS070 

ALMGS072 

(1.1.2) ), ALMGS074 

(1.1.3) ), ALMGS076 

(1.1.4) ), ALMGS078 

(1.1.5) ), ALMGS080 

(1.1.6) ), ALMGS082 

(1.1.7) ) ALMGS084 

ENTS ALMGS086 

ALMGS088 
EPOCH 1960 . ALMGS090 
EPOCH 1960. ALMGS092 
EPOCH 1960 .A LMGS094 
EPOCH I960. ALMGS096 
EPOCH 1965 .ALMGS098 
EPOCH 1965. ALMGS100 

ALTITUDE ALMGS102 
EPOCH 1970. ALMGS104 

48,99,3*0, ALMGS106 

0,-19498, ALMGS108 


2043,15853,12904,5026,2313,45,56,-88 ,74, 3*0 ,-4310 , 2308 , -1300 , 87 1 2ALMGS 1 10 
t-3940, -312,-2417, -75, -138, -156, 3*0, 15 20,-2684,29,-2505, 2714, ALMGS1 12 
-1573,-12,-244,-33,114,3*0,86, 12 12, -1160,-1 104, 799,-652, 5, -15,71, ALMGS 114 
111,3*0,-119,1028,609,-272,-124,-116,-1091,141,-56,10,3*0,-540, ALMGS 1 1 6 
-244,-91,22,276,-211,-201,58,117,4*0,69,-122,58,-170,26,236,-25, ALMGS 118 
-160,64,16,3*0,-220,156,51,-35,-18,96, 121,2,-25,15,42*0 / ALMGS1 20 
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,121, -1003, ALMGS122 
293,-924, ALMGS124 
*0,304, 288, ALMGS126 
ALMGS128 
ALMGS130 
5,104,-29, ALHGS132 
, —9 , ALMGS134 

,2*0,-4254, ALMGS136 
03,-2743, ALMGS138 
1 148,-1089, ALMGS140 
-133, -1089, ALMGS142 
, 16,-2,2*0,ALMGS144 
,-30,-19, ALMGS146 
20»11,28*0/ALMGS148 
1,2*0, -371, ALMGS150 
456,231, ALMGS152 
-4, 235,-90, ALMGS154 
,2*0,224, ALMGS156 
3,84, 23»-17ALMGS 158 
-50,-21,3, ALMGS160 
75,-46,31, ALMGS162 
ALMGS164 
43, 114,-18,ALMGS166 
,39,-27, 20, ALMGS168 
97,15,-11, ALMGS170 
,-11,15, ALMGS172 
,6, 1,2*0, ALMGS174r 
4, 19, -9, 4, ALMGS176 
ALMGS178 
4,110,-26, ALMGS180 
6, -20, -18, ALMGS182 

1.0, -4453, ALMGS184 
1354, -2667, ALMGS186 
-1287,-1151ALMGS188 
9,-43 ,-916,ALMGS190 
79,87,17, ALMGS192 

0, -204, 144, ALMGS194 

5,43,80, ALMGS196 

i,14*0/ ALMGS198 

1. 3.0, -466, ALMGS200 
0»-214,-441ALMGS202 
46,396,70, ALMGS204 

1, -20, 1 ,0,5,ALMGS206 
58,-80,50, ALMGS208 
132,523,81, ALMGS210 
03,-19,0, ALMGS2 1 2 
!6,-30,-3, ALMGS214 
17,22,-155, ALMGS216 

ALMGS2 1 8 
ALMGS220 
1,96,-17, ALMGS222 
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H-9, -4, 20, -11 ,1,9, -2, -2, 3, 4*0,- 
I -3, -4, -3, -3, -5, 56*0/ 

DATA GTT5 71,168*0/ 


UHlfl VJ I 1 I UtlOO* u/ 

DATA SHM I T ( 1 , 1 ) / 0.0 /♦ TMOLO / 0.0 /, MOD 
***** SUBSCRIPTED DO-LOOP VERSION BEGINS HERE 
DIMENSION CONST ( 13,13) , FN ( 13) ,FM( 13) 
DIMENSION P( 13,13)»DP( 13,13) ,SP( 13 ) , CP ( 13) 
DATA P( 1 ,1 ) ,CP( 1 ) ,DP( 1,1 ) ,SP( 1 ) / 2*1. ,2*0. 
** BEGIN PROGRAM 
1F( SHM IT (1,1).E0.-1.) GO TO 8 


DATA P( 1 ,1 ) ,CP( 1 ) ,DP( 1,1 ) ,SP( 1 ) / 
***** BEGIN PROGRAM 

1 F C SHM IT (1,1).E0.-1.) GO TO 8 

***** INITIALIZE * ONCE ONLY, FIRST 
SHM I T ( 1,1 )=-l. 

DO 18 N-1,13 

FN(N)=N 

DO 18 M=1 , 13 


r 1 ME SUBROUTINE IS CALLED 


ALMGS326 

ALMGS328 

ALMGS330 

ALMGS332 

ALMGS334 

ALMGS336 

ALMGS338 

ALMGS340 

ALMGS342 

ALMGS344 

ALMGS346 

ALMGS348 

ALMGS350 

ALMGS352 

ALMGS354 

ALMGS356 

ALMGS358 

ALMGS360 


L-3 



\ 


FM( M ) =M— 1 

18 CONST ( N, M ) = FLOAT ( ( N-2 } ** 2 — ( M— 1 )**2 ) / < ( 2*N-3 ) * ( 2*N-5 ) ) 

DO 2 N=2 *13 

SHMIT(N,1) = ( 2*N— 3 ) * SHMIT(N-1,1) / (N-l) 

J J = 2 

DO 2 M=2 * N 

SHM IT ( N »M ) = SHM I T ( N » M— 1 ) * SORT ( FLOAT (< N-M+ 1 ) *JJ >/( N+M-2 ) ) 
SHMIT(M-1,N)=SHMIT(N,M) 

2 JJ = 1 

DO 7 K=1 , 7 
F1 = L6( 1 ,1 »K ) 

F2=LGT ( 1 » 1 , K ) 

F3=LGTT ( 1 » 1 »K ) 

NMAX=NMX ( K ) 

L = 0 

DO 3 1 = 1 » 3 

3 ISUM( K, I ) = 0 
DO 4 N=1,NMAX 
DO 4 M= 1 »NMAX 
L = L+l 

ISUM(K,1)=ISUM(K,1 )+L*LG(N,M,K> 

I SUM( K, 2 ) = I SUM( K , 2 ) +L*LGT ( N » M , K ) 

4 I SUM( K»3 > = ISUM( K »3 )+L*LGTT ( N »M»K ) 

DO 6 1=1,3 

IF( I SUM ( K , I ) »EO.LSUM(K,I ) ) GO TO 6 
C ***** ERROR IN DATA CARDS - NOTE WRITE AND STOP STATEMENTS 
PRINT 5, K, I ,LSUM( K, I ) , ISUMI K, I ) 

5 FORMAT ( // /29H DATA WRONG IN ALLMAG — MODEL ,I2,3X,2HI=, II, 3X, 
A17HPRECALCULATED SUM , 1 10 , 3X , 1 7HTH I S MACHINE GETS, 110) 

STOP 

6 CONTINUE 

DO 7 N=1 , NMAX 
DO 7 M= 1 ,NMAX 

GG( N,M,K)=LG( N, M, K ) *SHM I T ( N,M ) / F 1 
GGT<N,M,K)=LGT( N,M,K )*SHMIT(N,M)/F2 

7 GGTT(N,M,K)=LGTT(N,M,K)*SHMIT(N,M>/F3 

8 IF( (MODEL. EO.MODOLD). AND. (TM.EQ.TMOLD) ) GO TO 11 
C ***** note WRITE STATEMENT - NEW MODEL OR NEW TIME 

PRINT 9, MODEL, (LABEL( I , MODEL), 1 = 1, 4), TM 

9 FORMAT ( ' 0 MODEL USED IS NUMBER ' , I 2 , 2X , 4A8 , • FOR TM =',F9.3/) 
IF( MODEL.LT .1 .OR. MODEL .GT.7 ) STOP 

MODOLD=MODEL 
TMOLD=TM 
NMAX=NMX( MODEL ) 

T=TM-TO ( MODEL ) 

DO 10 N= 1 , NMAX 
DO 10 M=1 ,NMAX 

10 G(N,M)=GG(N,M,MODEL)+T*(GGT(N,M,MODEL)+GGTT(N,M,MODEL)*T) - 
C ***** CALCULATION USUALLY BEGINS HERE 

11 SP( 2 )=SPH 
CP(2)=CPH 

DO 12 M= 3 , NMAX 

SP(M)=SP(2 )*CP( M-l )+CP( 2 )*SP(M-1 ) 

12 CP(M)=CP(2)*CP( M- 1 )— SP(2)*SP(M-1 ) 

A0R=6371.2/RKM 

AR=A0R**2 
BR=0.0 
BT=0.0 
BP=0 .0 

DO 17 N=2 ,NMAX 

L-4 


ALMGS362 

ALMGS364 

ALMGS366 

ALMGS368 

ALMGS370 

ALMGS372 

ALMGS374 

ALMGS376 

ALMGS378 

ALMGS380 

ALMGS382 

ALMGS384 

ALMGS386 

ALMGS388 

A L MG S3 90 

ALMGS392 

ALMGS394 

ALMGS396 

ALMGS398 

ALMGS400 

ALMGS402 

ALMGS404 

ALMGS406 

ALMGS408 

ALMGS410 

ALMGS41 2 

ALMGS414 

ALMGS416 

ALMGS418 

ALMGS420 

ALMG5422 

ALMGS424 

ALMGS426 

ALMGS428 

ALMGS430 

ALMGS432 

ALMGS434 

ALMGS436 

ALMGS438 

ALMGS440 

ALMGS442 

ALMGS444 

ALMGS446 

ALMGS448 

ALMGS450 

ALMGS452 

ALMGS454 

ALMGS456 

ALMGS458 

ALMGS460 

ALMGS462 

ALMGS464 

ALMGS466 

ALMGS468 

ALMGS470 

ALMGS472 

ALMGS474 

ALMGS476 

ALMGS478 

ALMGS480 



AR~AOR*AR 
DO 17 M=1 *N 
IF(M.EQ.N) GO TO 13 

P(NtMJ=CT*P(N— 1*M ) -CONST ( N *M ) *P ( N-2 *M ) 

DP(N f M)=CT«DP( N-1,M ) -ST*P< N-l f M ) -CONST ( N , M ) *DP ( N-2 , M ) 
GO TO 14 

13 P<NtN)=ST*P(N-ltN-l ) 

DPCNtN)~ST*DP<N-l,N-l )+CT*P ( N-l t N-1 ) 

14 PAR=P(N,M)*AR 

I F ( M. EQ. 1 ) GO TO 15 
TEMP-G( NfM)*CP(M) +G (M-lfN)*SP(M) 
BP*=BP-(G(NtM>*SP(M)~G(M-l, N)*CP(M) )*FM(M)*PAR 
GO TO 16 

15 TEMP = G(N,M) 

16 BR-BR-T EMP*FN(N)*PAR 

17 BT*BT+TEMP*DP<N*M)*AR 


ALMGS482 
ALMGS484 
ALMGS486 
ALMGS488 
ALMGS490 
ALMGS492 
ALMGS494 
ALMGS496 
ALMGS498 
ALMGS500 
ALMGS502 
ALMGS504 
ALMGS506 
ALMGS508 
ALMGS510 
ALMGS51 2 


1 BR s 

BR / 

100000. 

ALMGS514 

BT = 

BT / 

100000. 

ALMGS516 

BP = 

BP / 

ST / 100000. 

ALMGS518 

B * 

SQRT(BR*BR+8T*BT+BP*BP ) 

ALMGS520 


FOR DOUBLE PRECISION REPLACE PRECEDING CARD WITH FOLLOWING CARD ALMGS522 

B * DSQRT(BR*BR+BT*BT+BP*BP ) ALMGS524 

THIS AND THE IMPLICIT REAL*8 CARD ARE THE ONLY TWO CHANGES NEEDED* ALMGS526 
RETURN ALMGS528 

END ALMGS530 


0265 CARDS 
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SUBROUTINE DEK.MAG< MODEL , TM , RKM , ST , CT , SPH , CPH , BR , BT , BP , B ) 
it THIS VERSION READS IN COEFFICIENTS AS DATA CARDS. SEE BELOW. 

‘ GEOCENTRIC VERSION OF GEOMAGNETIC FIELD ROUTINE 
‘ SINGLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC, 029 PUNCH) 

‘ SHORT DECK, USES SUBSCRIPTED VARIABLES AND DO LOOPS 
‘ EXECUTION TIME PER CALL 3 TIMES GREATER THAN LONG DECK 
‘ PROGRAM DESIGNED AND TESTED BY E G ST ASS I NOPOULOS AND G D MEAD, 
‘ CODE 641, NASA GODDARD SPACE FLT CTR, GREENBELT, MD 20771 
INPUT: MODEL CHOICE .OF 7 MODELS 

'** TM TIME IN YEARS FOR DESIRED FIELD 

•** RKM GEOCENTRIC DISTANCE IN KILOMETERS 

'** ST ,CT SIN £ COS OF GEOCENTRIC CGLATITUDE 

‘** SPH, CPH SIN £ COS OF EAST LONGITUDE 


OUTPUT: BR ,BT ,BP GEOCENTRIC FIELD COMPONENTS 


GAUSS 


FIELD MAGNITUDE 


GAUSS 


NOTE: FOR GREATEST EFFICIENCY, COMPLETE ALL CALCULATIONS 


ONE MODEL AND ONE TIME BEFORE CHANGING MODELS OR TIME. 
***** FOR DOUBLE PRECISION ADD THE FOLLOWING CARD 
IMPLICIT REAL*8(A-H,0-Z) 

***** SEE END OF DECK FOR ONE MORE CHANGE 
DIMENSION T0( 7) ,NMX( 7) ,G( 13,13) 

REAL*8 LABEL (8,7) 

REAL*4 GG( 13, 13,7 ),GGT( 13, 13,7 ),GGTT( 13, 13,7), SHMIT( 13,13) 

DATA SHM I T ( 1 , 1 ) / 0.0 /, TMOLD / 0.0 /, MODOLD /Of 
DIMENSION CONST ( 13,13),FN(13),FM(13> 

DIMENSION P( 13,13 ) ,DP( 13,13 ) ,SP( 13) ,CP( 13) 

DATA P( 1,1 ),CP( 1) ,DP( 1,1), SP( 1) / 2*1 * , 2*0. / 

***** BEGIN PROGRAM 

I F ( SHM I T ( ltl).EQ.-l. ) GO TO 8 

***** INITIALIZE * ONCE ONLY, FIRST TIME SUBROUTINE IS CALLED 
SHM I T ( 1,1 >=-l. 

DO 18 N=1 ,13 
FN(N)=N 
DO 18 M=1 ,13 
FM(M.)=M-1 

CONST ( N ,M ) = FLOAT( ( N-2 ) **2-< M-l ) **2 ) / ( ( 2*N-3 ) * ( 2*N-5 ) ) 

DO 18 K= 1 , 7 
GG(N,M,K) = 0. 

GGT ( N, M, K ) = 0. 

18 GGTT ( N ,M ,K ) = 0. 

DO 2 N=2 , 1 3 

SHM I T ( N , 1 ) = ( 2*N— 3 ) * SHMIT(N-1,1) / (N-l) 

J J = 2 

DO 2 M=2 ,N 

SHMI T ( N, M ) = SHM I T ( N , M— 1 ) * SORT ( FLOAT (( N-M + 1 ) *JJ )/( N+M-2 ) ) 

SHMIT (M— 1,N)=SHMIT (N,M) 

2 JJ = 1 

COEFFICIENTS ARE READ IN WHEN |>£KMAG IS CALLED THE FIRST TIME 
SET UP INPUT DECK AS FOLLOWS: 

CARD 1. NMODLS ,NPR I NT (215). NMODLS = NO. OF MODLS TO BE READ IN 
DECK A: FIRST MODEL IN STANDARD CAIN FORMAT. 

DECK B: SECOND MODEL, ETC. 

FIRST CARD IN EACH DECK CONTAINS EPOCH (COLS 4-9) AND LABEL (10-73), 
LAST CARD IN EACH MODEL DECK CONTAINS ZERO OR BLANK IN COLS l-^. 
READ (5,4) NMODLS, NPRINT 
4 FORMAT (215) 

I F ( NMODLS. GT . 7 ) NMODLS = 7 
DO 3 K=l, NMODLS 

READ ( 5 , 20 ) TZ,(LABEL( I, K), 1=1, 8) 

20 FORMAT( 3X,F6.1 ,8A8) 


IS CALLED 


FIRST TIME 


MODLS TO BE 
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WRITE(6 t 21) K t (LABELIItK), 1=1,8) 

21 FORMAT( '0 MODEL NUMBER ■ 1 3 T 5X t 8 A8 ) 

MAXN=0 

5 READ (5t6) N »M t GNM tHNM t GTNM t HTNM t GTTNM t HTTNM 

6 FORMAT ( 2 1 3 1 6F1 1 * A ) 

IF (N.LE.O) GOT 07 

IF(NPRINT*EQ*1 ) WR I TE ( 6 , 6 ) N t M, GNM » HNM, GTNM, HTNM t GTTNM, HTTNM 
MAXN= (MAXO(NtMAXN) ) 

GG( Nt M t K ) = GNM 
GGT(NtMtK) = GTNM 
GGTT ( N*M t K ) = GTTNM 
IF (M.EQ.l) G0T05 
GG( M-l *N t K ) = HNM 
GGT ( M-l *N t K ) = HTNM 
GGTT(M-ltNtK) = HTTNM 
GO TO 5 

7 CONTINUE - 

IF(NPRINT.EQ.O) WRITE<6»22) GG(2tltK> 

22 FORMAT ( • G<2tl > ='F9*1 ) 

NMX(K) = MAXN 

TO ( K ) = TZ 
DO 3 N= 1 tMAXN 
DO 3 M= 1 t MAXN 

GG( N ?M t K ) = GG(NtMtK) * SHMIT(NtM) 

GGT < Nt Mt K ) = GGT(NtMtK) * SHM I T ( N, M ) 

3 GGTT(NtMtK) = GGTT(NtMtK) * SHM I T ( N t M ) 

8 I F ( ( MODEL •E0*M0D0LD)*AND*(TM*EQ*TM0LD) ) GO TO 11 
C ***** NOTE WRITE STATEMENT - NEW MODEL OR NEW TIME 

PRINT 9t MODELt <LABEL< I , MODEL ) t I =1 t 8 ) t TM 

9 FORMAT ( 1 0 MODEL USED IS NUMBER • 1 1 2 1 2X 1 8 A 8 t 1 FOR TM =*tF9.3/> 
IF( MODEL * LToUOR, MODEL. GT.NMODLS) STOP 

MODOLD = MODE L 

TMOLD=TM 

NMAX = NMX( MODEL ) 

T=TM-T0< MODEL > 

DO 10 N = 1 t NMAX 
DO 10 M= 1 1 NMAX 

10 G(N t M)=GG( NtMtMODEL )+T*( GGT (NtMtMODEL )+GGTT( N r Mt MODEL )«T) 

C ***** CALCULATION USUALLY BEGINS HERE 

11 SP { 2 ) =S PH 
CP( 2)=CPH 

DO 12 M = 3 t NMAX 

SP(M)=SP(2)*CP< M- 1 ) +CP ( 2 ) *SP < M- 1 ) 

12 CP<M)=CP<2)*CP(M-1 > — S P ( 2 )*SP(M — 1 ) 

AQR=6371*2/RKM 

AR=A0R**2 
BR = 0 *0 

BT = 0 *0 - 

BP=0* 0 

DO 17 N = 2 t NMAX 
AR=AOR*AR 
DO 17 M = 1 1 N 
IF(M.EO.N) GO TO 13 

P(N,M )=CT*P( N-l ,M ) -CONST (NtM)*P(N-2tM) 

DP(NtM) = CT^DP< N-l tM)-ST*P(N-ltM)-CONST(N,M )*DP(N-2,M) 

GO TO 1 A 

13 P(NtN)=ST*P{ N-l t N-l ) 

DP(NtN)=ST*DP<N-l f N-l )+CT*P<N-ltN-l ) 

1 A P AR= P ( N t M ) *AR 

IF(M.EQ.l) GO TO 15 



DEKMG122 

DEKMG12A 

DEKMG126 

DEKMG128 

DEKMG130 

DEKMG132 

DEKMG13A 

DEKMG136 

DEKMG138 

DEKMG1A0 

DEKMG1A2 

DEKMG1AA 

DEKMG1A6 

DEKMG1A8 

DEKMG150 

DEKMG152 

DEKMG15A 

DEKMG156 

DEKMG158 

DEKMG160 

DEKMG162 

DEKMG16A 

DEKMG166 

DEKMG1 68 

DEKMG1 70 

DEKMG172 

DEKMG17A 

DEKMG176 

DEKMG178 

DEKMG1 80 

DEKMG182 

DEKMG18A 

DEKMG186 

DEKMG1 88 

DEKMG190 

DEKMG192 

DEKMG19A 

DEKMG1 96 

DEKMG198 

DEKMG200 

DEKMG202 

DEKMG20A 

DEKMG206 

DEKMG208 

DEKMG210 

DEKMG2 12 

DEKMG21A 

DEKMG2 16 

DEKMG21 8 

DEKMG220 

DEKMG222 

DEKMG22A 

DEKMG226 

DEKMG228 

DEKMG230 

DEKMG232 

DEKMG23A 

DEKMG236 

DEKMG238 

DEKMG2A0 


I 


TEMP=G(N,M)*CP(M)+G(M-l t N)*SP(M) 
BP=BP-(G(N,M)*SP(M)-G(M-1,N)*CP(M) )*FM(M)*PAR 
GO TO 16 

15 TEMP = G( N *M ) 

16 BR=BR-TEMP*FN(N>*PAR 

17 BT=BT+TEMP*DP(N,M)*AR 
1 BR = BR / lOOOOO. 

BT = BT / 100000. 

BP = BP / ST / 100000. 

B = SQRT(BR*BR+BT*BT+BP*BP ) 

FOR DOUBLE PRECISION REPLACE PRECEDING CARD WITH FOLLOWING CARD 
B = DS0RT( BR*BR+BT*BT+BP*BP ) 

THIS AND THE IMPLICIT REAL*8 CARD ARE THE ONLY TWO CHANGES NEEDED 
RETURN 
END 


DEKMG242 

DEKMG244 

DEKMG246 

DEKMG248 

DEKMG250 

DEKMG252 

DEKMG254 

DEKMG256 

DEKMG258 

DEKMG260 

DEKMG262 

DEKMG264 

DEKMG266 

DEKMG268 

DEKMG270 


0135 CARDS 


SUBROUTINE ONEMAGI TM , RKM , ST , CT , SPH,CPH t BR , BT , BP» B ) 0NEMG002 

**** THIS VERSION CONTAINS ONE MODEL ONLY - IGRF 1965.0 0NEMG004 

**** GEOCENTRIC VERSION OF GEOMAGNETIC FIELD ROUTINE 0NEMG006 

**** SINGLE PRECISION DECK FOR IBM 360 MACHINES ( EBCDIC* 029 PUNCH) 0NEMG008 

**** SHORT DECK , USES SUBSCRIPTED VARIABLES AND 00 LOOPS ONEMGOIO 

**** PROGRAM DESIGNED AND TESTED BY G D MEAD, CODE 641, GSFC 0NEMG012 

***** INPUT: TM TIME IN YEARS FOR DESIRED FIELD 0NEMG014 

***** RKM GEOCENTRIC DISTANCE IN KILOMETERS 0NEMG016 

***** ST ,CT SIN £ COS OF GEOCENTRIC COLATITUDE 0NEMG018 

***** SPH»CPH SIN £ COS OF EAST LONGITUDE 0NEMG020 

***** OUTPUT: BR,BT,BP GEOCENTRIC FIELD COMPONENTS IN GAUSS 0NEMG022 

***** B FIELD MAGNITUDE IN GAUSS 0NEMG024 

***** FOR DOUBLE PRECISION ADD THE FOLLOWING CARD 0NEMG026 

IMPLICIT REAL*8(A-H, 0— Z ) 0NEMG028 

***** SEE END OF DECK FOR ONE MORE CHANGE 0NEMG030 

DIMENSION LG( 13, 13 ) , LGT ( 13 , 13 ) , G ( 1 3, 13 ) , GG ( 1 3, 1 3 ) , GGT < 13, 13), 0NEMG032 

1 SHMI T ( 13,13) 0NEMG034 

EQUIVALENCE I LGI 1 , 1 ) ,GG( 1 , 1 ) ) , I LGT 1 1 , 1 ) , GGT I 1 , 1 ) > 0NEMG036 

DATA SHMIT ( 1 ,1 )/0./»TM0LD/0./»TZER0/1965./»NMAX/9/ 0NEMG038 

DATA LG / 1, -30339,-1654,1297,958,-223,47,71,10,4*0,5758,-2123, 0NEMG040 
A 2994,-2036,805,357,60,-54,9,4*0,-2006,130,1567,1289,492,246,4,0, 0NEMG042 
B -3, 4*0, -403, 242 ,-176, 843, -392, -26, -229, 12,-12,4*0, 149,-280,8,-2650NEMG044 
C ,256,-161 ,3, -2 5, -4, 4*0, 16,125,-123,-107,77,-51,-4,-9,7,4*0,-14, 0NEMG046 
D 106,68,-32,-10,-13,-112,13,-5,4*0,-57,-27,-8,9,23,-19,-17,-2, 12, 0NEMG048 
E 4*0,3,-13,5,-17,4,22,-3,-16,6,56*0/ 0NEMG050 

DATA LGT / 10, 1 53,-244, 2 , -7 , 19, -1 , -5 , 1 , 4*0, -2 3, 87 , 3,- 108 , 2 , 1 1 , -3, 0NEMG052 
F -3,4,4*0,-118,-167,-16,7,-30,29,11,-7,6,4*0,42,7,-77,-38,-1,6, 19,0NEMG054 
G -5,5*0,-1,16,29,-42,-21,0,-4,3,5*0,23, 17, -24, 8, -3, 13, -4, 0,-1, 4*0, 0NEMG056 
H— 9 ,— 4 ,20 ,— 1 1 ,1,9,— 2 ,— 2 ,3,4*0,— 11,3,4,2,4,2,3,— 6, —3 ,4*0,1,— 2,— 3,-2, 0NEMG058 
I -3, -4, -3, -3, -5, 56*0/ QNEMG060 

DIMENSION CONST! 13,13), FN( 13), FM( 13) 0NEMG062 

DIMENSION P( 13,13) , DPI 13,13) ,SP( 13 ) ,CP( 13) 0NEMG064 

DATA PI 1 ,1 ) ,CP( 1 ) ,DP( 1,1 ) ,SP( 1) / 2*1. ,2*0. / 0NEMG066 

IFISHMITI 1,1).EQ.-1. ) GO TO 8 0NEMG068 

***** INITIALIZE * ONCE ONLY, FIRST TIME SUBROUTINE IS CALLED 0NEMG070 

SHMITI 1, 1 )=-l. 0NEMG072 

DO 18 N=1 , 1 3 0NEMG074 

FNI N ) =N 0NEMG076 

DO 18 M=1 ,13 0NEMG078 

FMI M ) =M— 1 0NEMG080 

1 8 CONST I N ,M ) = FLOAT I I N-2 ) **2— I M— 1 ) **2 ) / II 2*N-3 ) * I 2*N-5 ) ) 0NEMG082 

DO 2 N=2» 13 0NEMG084 

SHMIT I N , 1 ) = (2*N— 3) * SHMITIN-1,1) / IN-1) 0NEMG086 

JJ=2 0NEMG088 

DO 2 M=2,N 0NEMG090 

SHMIT I N, M ) = SHM I T I N »M-1 ) * SORT I FLOAT II N-M+l ) * J J ) / I N+M-2 > ) 0NEMG092 

SHMIT ( M— 1 , N ) =SHMIT I N ,M ) 0NEMG094 

2 JJ = 1 ~ ' 0NEMG096 

FI = LGI 1 , 1 ) 0NEMG098 

F2 = LGT (1,1) 0NEMG100 

DO 7 N=1 ,NMAX 0NEMG102 

DO 7 M= 1 , NMAX 0NEMG104 

GGI N,M ) = LGI N ,M ) *SHM IT I N ,M ) /F 1 0NEMG106 

7 GGT I N, M ) = LGT I N, M ) *SHM IT(N,M)/F2 0NEMG108 

8 IFITM.EQ.TMOLD) GO TO 11 0NEMG110 

TMOLD=TM 0NEMG112 

T = TM - TZERO 0NEMG114 

DO 10 N= 1 , NMAX 0NEMG116 

DO 10 M=1,NMAX 0NEMG1 1 8 

10 G(N,M) = GGI N, M ) + T*GGT I N,M ) 0NEMG120 


FLOAT! (N-2)**2-(M-l)**2) / I I 2*N-3 ) * I 2*N-5 ) ) 


L-9 


c ***** CALCULATION USUALLY BEGINS HERE 

11 SP ( 2 ) = SPH 
CP ( 2 ) =CPH 

DO 12 M=3 1 NMAX 

SP(M)=SP(2)*CP( M— 1 ) +CP ( 2 ) *SP ( M-l ) 

12 CP<M)=CP< 2)*CP<M-1 )-SP(2)*SP(M-l ) 

AOR=6371.2/RKM 

AR=AOR**2 
BR-0 .0 
BT=0.0 
BP=0.0 

DO 17 N=2 tNMAX 
AR=AOR*AR 
DO 17 M= 1 * N 
IFCM.EQ.N) GO TO 13 

P(NtM)=CT*P<N-l,M)-C0NST(N t M)*P<N-2tM) 

DP(N»M)=CT*DP(N-1 t M)-ST*P( N-l t M) -CONST IN, M)*DP(N-2,M) 

GO TO 14 

13 P(N,N) = ST*P(N-l r N-l ) 

DP(N*N )=ST*DP( N— 1 ?N— 1 )+CT*P(N-l *N-1 ) 

1 4 PAR = P ( N »M ) *AR 

I F ( M. EQ. 1 ) GO TO 15 
TEMP=G(N,M)*CP(M)+G(M-1,N)*SP(M) 

BP-BP— (G(N*M)*SP(M)-G< M-l ,N)*CP(M) )*FM<M)*PAR 
GO TO 16 

15 TEMP = G(NtM) 

16 BR=BR-TEMP*FN<N)*PAR 

17 BT-BT+TEMP*DP(N*M)*AR 
BP = BP/ST/100000, 

BR = BR/ 100000, 

BT = BT/100000. 

B * SORT(BR*BR+BT*BT+BP*BP ) 

C FOR DOUBLE PRECISION REPLACE PRECEDING CARD WITH FOLLOWING CARD 
C B * DSGRT(BR*BR+BT*BT+BP*BP ) 

C THIS AND THE IMPLICIT REAL*8 CARD ARE THE ONLY TWO CHANGES NEEDED. 
RETURN 
END 


0NEMG122 
0NEMG1 24 
0NEMG1 26 
0NEMG1 28 
0NEMG130 
0NEMG132 
0NEMG1 34 
0NEMG1 36 
ONEMG138 
0NEMG140 
0NEMG142 
0NEMG144 
0NEMG146 
0NEMG148 
0NEMG150 
0NEMG152 
0NEMG1 54 
0NEMG156 
ONEMG158 
0NEMG1 60 
0NEMG1 62 
0NEMG164 
0NEMG166 
0NEMG168 
0NEMG170 
0NEMG1 72 
0NEMG1 74 
0NEMG1 76 
ONEMG178 
0NEMG1 80 
0NEMG1 82 
0NEMG1 84 
0NEMG186 
0NEMG1 88 
0NEMG190 
0NEMG192 
0NEMG1 94 

0097 CARDS 
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c 

**** STANDARD CALLING 

ROUTINE 

FOR GDALMG 


002 

c 

**** SINGLE PRECISION 

DECK FOR 

IBM 360 MACHINES 

(EBCDIC, 029 PUNCH! 

004 

c 

**** NEEDS SINGLE-PRECISION SUBROUTINES GDALMG 

AND ALLMAG 

006 


RE AL*4 TM( 2 ) / 1965. 

f 1970. / t 

GDLAT( 2)/-30. ,60./ 

,GLON ( 2 ) /-90. ,90./, 

008 


1 ALT/100./ 




010 


PRINT 1G 




012 


10 FORMAT ( '1 MODEL 

TIME 

LAT LONG 

ALT X Y 

014 


1 z 

F 

H DEC 

INC') 

. 016 

* 

DO 20 M0DEL=1 ,7 




018 


DO 20 1=1,2 




020 


DO 20 J=1 ,2 



» 

022 


00 20 K=l,2 




024 


CALL GDALMG(MODEL,TM< I ) ,GDLAT< J).GLON(K) , ALT , X , Y , Z , F, H, DEC , A I NC ) 026 

20 PRINT 30, MODEL, TM( I ),GDLAT( J ) , GLON ( K ) , ALT , X , Y, Z , F , H, DEC, AI NC 028 

* 30 FORMAT! I 5 , 3X ,4F8 .0, 5F 10. 5 ,2F 10 .2 ) 030 

STOP 032 

END 034 


0017 CARDS 
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SUBROUTINE GDALMG( MODEL , TM ,GDL AT , GDLQN, GDALT t X , Y , Z , F , H, DEC , A I NC ) 

♦ ♦♦♦ GEODETIC VERSION OF GEOMAGNETIC FIELD SUBROUTINE 

♦♦♦♦ SINGLE PRECISION DECK FOR IBM 360 MACHINES ( EBCDIC? 029 PUNCH) 

♦ ♦♦♦ NEEDS SINGLE-PRECISION SUBROUTINE ALLMAG 

♦♦♦♦ PROGRAM DESIGNED AND TESTED BY G D MEAD AND E G ST ASS I NOPOULOS , 

♦ ♦♦♦ CODE 641, NASA GODDARD SPACE FLT CTR, GREENBELT, MD 20771 

♦♦♦INPUT: MODEL = 1-7; CHOICE OF SEVEN MODELS - SEE ALLMAG 

TM = TIME IN YEARS FOR DESIRED FIELD < E.G. 1971*25) 

GDLAT = GEODETIC LATITUDE (DEGREES) 

GDLON = EAST LONGITUDE (DEGREES) 

GDALT = ALTITUDE ABOVE GEO I D { KMS ) 

♦♦♦OUTPUT: X , Y , Z - GEODETIC FIELD COMPONENTS (GAUSS) 

F * MAGNITUDE OF FIELD (GAUSS) 

H = HORIZONTAL INTENSITY (GAUSS) 

DEC, AINC = DECLINATION AND INCLINATION ANGLES (DEGREES) 

♦♦♦ NOTE: FOR GREATEST EFFICIENCY, COMPLETE ALL CALCULATIONS WITH ONE 
MODEL AND ONE TIME BEFORE CHANGING MODEL OR TIME. 

REFERENCE GEOID IS THAT ADOPTED BY IAU IN 1964 

DATA RAD, A, AB2,E2/ 57. 29578, 6378. 16, 1.0067397,. 0067397/ 

SINLAT = S IN ( GDLAT/RAD ) 

COSLAT = SORT ( l.-SINLAT*^2) 

I F ( MODEL • EQ. 6 ) GO TO 2 

1 SINBET * SINLAT / SQRT(SINLAT**2+( AB2^C0SLAT ) **2 ) 

:♦♦♦ BETA = GEOCENTRIC LATITUDE AT SURFACE OF GEOID 

COSBET = SQRT ( 1 .-S INBET^S INBET ) 

RGEOID * A / SQRT( 1 . +E2*S I NBET*S I NBET ) 

XKM = RGEOID ♦COSBET + GDALT *C0 SL AT 
YKM = RGEOID^SINBET + GDALT^S I NL AT 
RKM a SORT ( XKM^«2+YKM**2 ) 

ST = XKM/RKM 
CT = YKM/RKM 
GO TO 3 

2 RKM = 6371.2 + GDALT 
ST = COSLAT 

CT = SINLAT 

3 SP = S IN ( GDLON/RAD ) 

CP = COS (GDLON/RAD) 

CALL ALLMAG(MODEL,TM,RKM,ST,CT,SP,CP,BR,BT, Y,F ) 

SIND = ST*S I NL AT - CT^COSLAT 
COSD = CT*S I NLAT + ST*COSLAT 
X = -BT ♦COSD - BR^SIND 
Z = BT^SIND - BR^COSD 
H * SQRT ( X*X+Y*Y ) 

DEC = RAD ♦ ATAN2 ( Y, X ) 

AINC = RAD ♦ AT AN ( Z/H ) 

RETURN 

END 


GDALT 


GDALM002 

GDALM004 

GDALM006 

GDALM008 

GDALM010 

GDALMO 1 2 

GDALMO 14 

GDALMO 1 6 

GDALMO 18 

GDALM020 

GDALM022 

GDALM024 

GDALM026 

GDALM028 

GDALM030 

GDALM032 

GDALM034 

GDALM036 

GDALM038 

GDALM040 

GDALM042 

GDALM044 

GDALM046 

GDALM048 

GDALM050 

GDALM052 

GDALM054 

GDALM056 

GDALM058 

GDALM060 

GDA-LM0 62 

GDALM064 

GDALM066 

GDALM068 

GDALM070 

GDALM072 

GDALM074 

GDALM076 

GDALM078 

GDALM080 

GDALM082 

GDALM084 

GDALM086 

GDALM088 

GDALM090 

GDALM092 

GDALM094 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


.INPUT: 


* **** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 

*****OUTPUT 

***** 

***** 

***** 

***** 

***** 


MAIN DECK FOR LINTRA LINTR002 
MODIFIED FIELD LINE TRACING ROUTINE AS OF JAN 1971, LINTR004 
DESIGNED AND TESTED BY E G ST ASS I NOPOULOS AND G D MEAD, LINTR006 
CODE 641, NASA GODDARD SPACE FLT CTR, GREENBELT, MD 20771 LINTR008 


SINGLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC, 029 PUNCH ) L 1NTR010 
NEEDS SINGLE-PRECISION SUBROUTINES CONVRT, I TER AT, AND ALLMAG. L INTR01 2 
ACCEPTING INPUT COORDINATES OF STARTING POINTS IN EITHER A LINTR014 
GEOCENTRIC OR GEODETIC SYSTEM AND RETURNING COORDINATES LINTR016 
FOR CONJUGATE INTERSECT POSITIONS IN BOTH SYSTEMS. LINTR018 
INTERSECTS MAY BE OBTAINED AT ANY SPECIFIED ALTITUDE tEVEL. LINTR020 
DIRECTION OF TRACING INPUT-CONTROLLED THROUGH PARAMETER "D I R" .L I NTR022 
FOR SURFACE POINTS "DIR" SHOULD ALWAYS BE <+l). LINTR024 
FOR SPACE POINTS, <DIR>0) WILL TRACE THE FIELD LINE TO THE LINTR026 
OPPOSITE HEMISPHERE, (DIRCO) WILL TRACE THE FIELD LINE TOWARDSL INTRO 28 
THE SURFACE IN THE SAME HEMISPHERE. 


MODEL 

TM 

NPRINT 

ICOORD 

GDLAT ,GLON,ALT 
GCLAT ,GLON,RKM 
DS 
DIR 


CHOICE OF 7 MODELS (FROM ALLMAG) 

TIME IN YEARS FOR DESIRED FIELD 

PARAMETER CONTROLLING OUTPUT 

=0 NO RUNNING PRINTOUT; =1 PRINT STEPS 

REFERENCE SYSTEM OF INPUT COOROINTATES 

=1 GEODETIC; =2 GEOCENTRIC 

GEODETIC STARTING POINT COORD ( OEGR. , KM ) 

II 


LINTR030 

LINTR032 

LINTR034 

LINTR036 

LINTR038 

LINTR040 

LINTR042 


HALT 

LABEL 

PLAT , PLON , PRKM 
PGLAT t PL ON » P AL T 
ARC 
BMIN 

BM I NLT t BM I NLN i 
BMINR 


LINTR044 

GEOCENTRIC STARTING POINT COORD( » )LINTR046 

TRACING STEPSIZE IN KILOMETERS LINTR048 

PARAMETER CONTROLLING DIRECTION OF TR AC EL I NTR050 
+1. STARTS TRACING TOWARDS HIGHER ALT* LINTR052 
-1. STARTS TRACING TOWARDS LOWER ALT . LINTR054 

GEODETIC ALTITUDE OF CONJUGATE I NTER SECTL I NTR05 6 
NAME OF ORIGIN (STARTING POINT) LINTR058 

GEOCENTRIC COORD* OF CONJUGATE I NTER SECTL I NTR060 
GEODETIC COORD. OF CONJUGATE INTERSECT LINTR062 


ARCLENGTH OF FIELD LINE TRACED, IN KM 
MINIMUM FIELD STRENGTH ALONG LINE 
GEOCENTRIC COORD. OF BMIN POSITIONt IN 
DEGREES AND KILOMETERS 


LINTR064 
LINTR066 
LINTR068 
LINTROTO 
LINTR072 
LINTR074 
L INTR076 


COMMON /ITER/ L * R , DL AT t DLON , RP , DL ATP * DLONP , BR , BT , BP , B, ST, SGN, DS 
DATA RAD/57.2957795/, C 1 / .0067397 / , RA/6378.16/, MAXS/200/ 

10 FORMAT ( 6F10.6,2A4) 

11 FORMAT ( *0’ ,47X , 1 S T E P • , 7X , • L AT ■ , 5X , 1 LON 1 ,4X , * RKM » , 3X , • ALT • , 7X , * BR • ,LINTR078 

16X, 1 BT ' , 6X, ■ BP* , 7X, 1 B 1 , / ) LINTR080 

12 FORMAT( 20X , 'OLD COORDINATES FOR STEP# 1 , I 6 , * ** 1 ,2F9.3,2F8.0,4F9.5)LINTR082 

13 FORMAT ( 1 1 MODEL ’tie/,' TIME ',F9.2/t' PRINT 


FORMAT ( 1 1 MODEL 1 ,18/, 1 
14 FORMAT( 10X, 'GEOCENTRIC 
1SIZE/ARCLENGTH DIR 
27 X , * L AT LONG ALT' 


1 , I 8 / , * COORD 1 , 18,/) LINTR084 


COORDINATES GEODETIC COORDINATES 
HALT LABEL 1 ,/,12X, 'LAT LONG RKM 
,9X, •DS/ARC 1 ,/,HX, ' (DEGR) (DEGR) (KM) 


15 

16 

17 

18 
19 


36X , ■ (DEGR) (DEGR) ( KM ) • ,9X , • < KM ) 1 , / ) 

FORMATt ' ORIGIN ' , 2 F8 .2 , F8 . 1 , 1 X2F 8 . 2 , F 8 . 1 , F 1 1 . 0 , 9X , 2 F7. 0 , 5X , 2 A4 ) 
FORMAT ( 15, F10. 2,215) 

FORMAT ( * 0 I NTRSCT 1 , F6 .2 , F8 . 2 * F8 . I , IX * 2F8 . 2 1 F8 . 1 » FI 1 .0* //// ) 
FORMATt • CHECK INPUT: ALT= • , F 8 . 0 , 9X , • D I R= ■ , F3. 0 , 9X, , HALT=',F8 


STEPL I NTR086 
*, LINTR088 
*, LINTR090 
LINTR092 
LINTR094 
LINT R096 
LINTR098 
0/ / ) L I NTR 100 


FORMAT ( 1 OL I NE TRACING TERMINATED: ITERATION EXCEEDS 200 STE PS 1 , / / ) L I NTR 1 02 
PINTER(A1,A2,A3,A4,A5,A6,A7) = ( < A2- A3 ) * ( A7-A2 ) * ( A7-A3 ) *A4- ( A 1-A3 ) L I NTR 104 

1 *t A7-A1 )*( A7-A3 ) *A 5 + ( A 1-A2 ) * ( A7-A 1 ) * ( A7-A2 ) *A6 ) / ( ( A 1 -A 2 ) * ( A 1 -A3 ) LI NTR 106 

2 *< A2-A3) ) L I NTR 108 

READ(5,16,END = 6) MODEL ,TM ,NPR I NT , I COORD LINTR110 

WR I TE ( 6f 1 3 ) MODEL, TM^NPRINT, ICOORD LINTR112 

IF(NPRINT.EQ.O) WRITE (6,14) LINTR114 

IF( ICOORD. EO. 2) GO TO 2 LINTR116 

READ( 5,10,END = 6 ) GDL AT ,GLON , ALT ,DS ,D1R , HALT , LABEL 1 ,LABEL2 L1NTR118 

CALL CONVRT( 1, GDLAT, ALT, GCLAT, RKM) LINTR120 


L-13 


GO TO 3 

READC5, 10,END=6) GCL AT , GLON , RKM , DS , D I R * HALT , L ABEL 1 , LABEL 2 
CALL CONVRT( 2 ,GDLAT ,ALT,GCLAT,RKM ) 

I F ( NPR I NT • EQ. 1 ) WRITE(6,14) 

SINGML-.98*SIN( GCLAT/RAD ) + . 199*C0 S ( GCLAT/R AD ) *COS ( (GLON+69. 
DS = .06*RKM/(1.-SINGML**2) - 370, 

I F ( DS .GT *3 • E3 ) DS=3.E3 

IF< ( DIR.LT.O. ).AND.(DS.GT. (ALT-HALT )/20. ) ) DS= ( ALT-HAL T ) /20 
WRITE (6, 15 )GCL AT , GLON, RKM ,GDLAT, GLON, ALT , DS , D I R , HALT , LABEL 1 
IF( ALT+DIR*HALT.GE. 0.0) GO TO 7 
WR I TE ( 6 » 18 ) ALT, DIR, HALT 
GO TO ( 1 , 2 ) t ICOORD 
I F ( NPR I NT • EQ • 1 > WRITE(6,11) 

BM IN=5 • 0 
DL AT = GCLAT 
DLON=GLON 
R = RKM 
L = 0 

CT=S I N ( DLAT/RAD ) 

ST=SORT ( 1 • -CT*CT ) 

SP=SIN( DLON/RAO) 

CP=COS(DLON/RAD) 

CALL ALLMAG(MODEL,TM,R,ST ,CT ,SP,CP,BR,BT ,BP,B) 

I F ( L • EO. 0 ) SGN=S IGN( l.,BR*DIR) 

I F ( BMIN.LE.B) GO TO 5 

BM I N=B 

BM INLT=DLAT 

BM INLN=DLON 

BM I NR=R 

L==L+ 1 

IF(L.LT.MAXS) GO TO 8 
WR I T E ( 6, 1 9 ) 

GO TO (1,2) , ICOORD 
DL ATPP=DL ATP 
DLONPP=DLONP 
CALL ITERAT 
HPP=HP 

HP=RP-RA/SQRT( l.+Cl*CT*CT ) 

I F ( NPR I NT • EO • 1 ) WR I T E ( 6 , 1 2 ) L , DL ATP , DLONP , RP , HP , BR , BT , BP , B 
IF(R.GT.RP) GO TO 4 

H=R-R A / SORT ( 1 . +C 1 *S I N ( DL AT/R AD ) **2 ) 

IF(H.GT.HALT) GO TO 4 

PLAT=PINTER( HP P ,HP , H , DL AT PP , DL ATP , D L AT ,H ALT ) 

PLON=P INTER ( HPP,HP,H,DLONPP,DLONP,DLON,HALT ) 

PRKM=HA LT+RA/ SORT (. 1 . +C 1*S I N ( PL AT /R AD ) **2 ) 

ARC=DS*( ( L-5 ) + ( HP-HALT )/<HP-H) > 

IF( PLON.LT.-180.0) PL ON =PL ON +360.0 
IF(PLON.GT. 180.0) PL 0N = PLON— 360 • 0 
CALL CONVRT (2,PGDLAT ,PGALT ,PLAT ,PRKM) 

WRITE( 6, 17) PLAT, PLON , PRKM , PGDL AT , PLON , PGAL T , ARC 

GO TO (1,2) , ICOORD 

STOP 

END 


LINTR122 
LINTR124 
LINTR126 
LINTR128 
) /RAD ) LINTR130 
LINTR132 
LINTR134 
. LINTR136 

,LABEL2LINTR138 
LINTR140 
LINTR142 
LINTR144 
LINTR146 
LINTR148 
LINTR150 
LINTR152 
LINTR154 
LINTR156 
LINTR158 
LINTR160 
LINTR162 
LINTR164 
LINTR166 
LINTR168 
LINTR170 
LINTR172 
LINTR174 
LINTR176 
LINTR178 
LINTR180 
LINTR182 
LINTR184 
LINTR186 
LINTR188 
LINTR190 
LINTR192 
LINTR194 
LINTR196 
LINTR198 
LINTR200 
LINTR202 
LINTR204 
LINTR206 
LINTR208 
LINTR210 
LINTR21 2 
LINTR214 
LINTR216 
LINTR218 
LINTR220 
LINTR22 2 
LINTR224 
LINTR226 
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SUBROUTINE ITERAT 

SINGLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC, 029 PUNCH) 
FIELD LINE INTEGRATION PROGRAM USING A 4-POINT ADAMS FORMULA AFTER 
INITIALIZATION. FIRST 7 ITERATIONS ADVANCE POINT BY. 3*DS 


C*** 



Y ( 1-3 ) : R, 

c*** 

B,BR,BT, BP 

C*** 

ST 


SGN 



c*** 

DS 

C***OUTPUT : 

Y( 1-3) : R, 


YOLD ( 1-3) : 


O** INPUT: L STEP COUNT. SET L=1 FIRST TIME THRU; 

:* SET L=L+ 1 THEREAFTER. 

<* Y ( 1-3 ) : R ,DL AT , DLON GEOCENTR TRACING POINT COORD .( KM , D6GR ) 

* B,BR,BT,BP FIELD t COMPONENTS AT POINT Y(l-3) 

ST SINE OF GEOCENTRIC COLATITUDE 

:* SGN SGN=+1 : TRACES IN DIRECTION OF FIELD 

SGN=-l: TRACES OPPOSITE TO FIELD DIREC 
;* DS INTERGR.STEPSIZE (ARC INCREMENT) IN KM 

^OUTPUT: Y ( 1-3 ) : R ,DL AT ,DLON NEW IMPLEMENTED TRACING POINT COORD 

YOLD ( 1-3 ) : RP , DL AT P , DLONP OLD Yll-3), BEFORE IMPLEMENTATION 
COMMON /ITER/ L , Y , YOLD , BR , BT , BP , B , ST , SGN , DS 
DIMENSION Y(3),Y0LD(3),YP(3,4) 

DATA RAD /57. 2957795/ 

YP ( 1,4) = SGN* ( BR /B ) 

FAC=SGN*RAD/(6*Y( 1 ) ) 

YP( 2,4 )=-BT*FAC 
YP( 3,4)=BP*FAC/ST 
IF( L.GT.7) GO TO 9 
DO 8 1=1,3 

GO T 0 ( 1,2, 3, 4, 5, 6, 7) ,L 

1 D2 = DS / 2 • 

D6=DS/6. 

D12=DS/ 12 • 

D24=DS/24. 

YP ( 1,1) = YP( I ,4) 

YOLD( I ) = Y( I ) 

Y{ I ) “= Y0LD( I ) + DS* YP( 1,1) 

GO TO 8 

2 YP( 1,2) = YP ( 1,4) ‘ 

Y ( I ) = YOLD(I) + D2 * ( YP ( 1,2) + YP(I,1)) 

GO TO 8 

3 Y ( I ) = YOLD(I) + D6 * <2.*YP(I,4) + YP(I,2) + 3.*YP(I,1)) 

GO TO 8 

4 YP( 1,2) = Y P ( 1,4) 

YOLD ( I) = Y ( I ) 

Y ( I ) = YOLD(I) + D2 * (3.*YP(I,2) - YP ( I , 1 ) ) 

GO TO 8 

5 V(I) = YOLD ( I ) + D12 * <5.*YP(I,4) + 8.*YP(I,2) - YP(I,1)) 

GO TO 8 

6 YP( 1,3) = YP( 1,4) 

YOLD ( I ) = Y ( I ) 

Y(I) = YOLD(I) + D12 * ( 23 • * YP ( 1,3) - 16.*YP(I,2) + 5#*YP(I f l)) 

GO TO 8 

7 Y( I ) = YOL D ( I ) +D24* ( 9. * YP ( I ,4)+19.*YP( I , 3 )-5.*YP ( 1,2) +YP (1,1)) 

8 CONTINUE 
RETURN 

9 DO 10 1=1,3 
YOLD ( I ) = Y( I ) 

Y ( I ) = YOLD ( I ) +D24* ( 55.*YP( I ,4)-59.*YP( I ,3 )+37.*YP( I ,2)-9.*YP( 1,1)) 
DO 10 J= 1 , 3 
10 YP ( I , J ) = YP ( I ,J+1 ) 

RETURN 
END 


- YP( 1,1) ) 


I TRAT002 

I TRAT004 

I TRAT006 

ITRAT008 

ITRAT010 

I TRAT01 2 

1 TR AT 014 

ITRAT016 

ITRAT018 

ITRAT020 

ITRAT022 

ITRAT024 

1TRAT026 

I TRAT028 . 

1TRAT030 

I TRAT032 

1TRAT034 

ITRAT036 

I TR AT 038 

ITRAT040 

ITRAT042 

I TRAT044 

ITRAT046 

I TRAT048 

ITRAT050 

I TRAT052 

1 TRAT054 

ITRAT056 

I TRAT058 

ITRAT060 

ITRAT062 

I TRAT064 

ITRAT066 

ITRAT068 

IT RAT 070 

ITRAT072 

ITRAT074 

ITRAT076 

I TR AT 078 

I TRAT080 

ITRAT082 

ITRAT084 

I TR AT 086 

I TRAT088 

1 TRAT090 

ITRAT092 

I T RATO 9 4 

ITRAT096 

I TRAT098 

ITRATIOO 

ITRAT102 

ITRAT104 

ITRAT106 

ITRAT108 

ITRAT1I0 

ITRAT112 

I TR AT 114 
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SUBROUTINE CONVRTI I , GDLAT , ALT , GCLAT , RKM ) 

C SINGLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC* 029 PUNCH) 

C CONVERTS SPACE POINT FROM GEODETIC TO GEOCENTRIC OR VICE VERSA 
C REFERENCE GEO I D IS THAT ADOPTED BY IAU IN 1964 
C A=6378.16, 8=6356.7746, F=l/298.25 
C I = I GEODETIC TO GEOCENTRIC 
C 1=2 GEOCENTRIC TO GEODETIC 
C GDLAT = GEODETIC LATITUDE IN DEGREES 
C ALT = ALTITUDE ABOVE GEOID IN KILOMETERS 
C GCLAT = GEOCENTRIC LATITUDE IN DEGREES 
C RKM = GEOCENTRIC DISTANCE IN KILOMETERS 

DATA A,RAD,AB2,EP2/6378. 16,57.29578, 1.0067397,. 0067397/ 

I F ( I.EQ.2) GO TO 2 

1 SINLAT = S IN( GDLAT /RAD ) 

COSLAT = SQRT( 1 ,-S INLAT**2 ) 

COSTH = SINLAT / SQRT< ( AB2*C0SLAT )*#2+SINLAT**2 ) 

SINTH = SORT ( 1 .-C0STH4*2 ) 

RGEOID = A / SQRT< 1 .+EP2*C0STH**2 ) 

X = RGEOID*SINTH + ALT *COSL AT 
Y = RGEO I D*COSTH + ALT = S I NL AT 
RKM = SORT ( X*X+Y*Y ) 

GCLAT = RAD * ATAN ( Y/X ) 

RETURN. 

2 RER=RKM/A 

C SEE ASTRON. J. VOL. 66 P.15, 1961, FOR FORMULAS BELOW. 

A2=< ( - 1.4127348E-8/RER+ .9433913 IE-8 )/RER+.33523288E-2)/RER 
A4= ( ( ( -1 .2 545063 E—10/RER+ • 1 1 76 0996 E -9 ) /RER+. 1 1 238084E-4 ) /RER 
1 -.28.14244E-5 ) / RER 

A6= ( ( 54.939685E-9/RER-28.301730E-9) /RER+3.5435979E— 9 ) /RER 
A8=( ( ( 320.00/RER-252.00 ) /RER+64.00 ) / RER- 5. 00 )/RER* c .98008304E-12 
SCL = S I N ( GCLAT/RAD ) 

CCL = SQRT< l.-SCL*SCL ) 

S2CL=2.*SCL*CCL 

C2CL=2.*CCL*CCL-1.0 

S4CL=2 »*S2CL*C2CL 

C4CL=2.*C2CL*C2CL-1.0 

S8CL=2.*S4CL*C4CL 

S6CL=S2CL*C4CL+C2CL*S4CL 

DLTCL=S2CL*A2+S4CL*A4+S6CL*A6+S8CL*A8 

GDLAT =DLTCL*RAD+ GCLAT 

ALT = RKM - A/SORT ( 1 .+EP2*SCL*SCL > 

RETURN 

END 

Following are Input data cards for the two LINTRA test runs: 

1 1965.0 1 1 


45.00 

-90.00 

o 

• 

o 

300.0 

1.0 

0 . 

TEST 

1 

1 

1965.0 0 

1 






75.00 

0.00 

0.0 

300.0 

1.0 

0 . 

TEST 

RUN 

45.00 

-90.00 

0.0 

300.0 

1.0 

0. 

TEST 

RUN 

-45.00 

90.00 

0.0 

300.0 

1.0 

0. 

TEST 

RUN 

-75.00 

0.00 

0.0 

300.0 

1.0 

0. 

TEST 

RUN 


CNVRT002 

CNVRT004 

CNVRT006 

CNVRT008 

CNVRT010 

CNVRTO 1 2 

CNVRT014 

CNVRTO 1 6 

CNVRT018 

CNVRT020 

CNVRT022 

CNVRT024 

CNVRT026 

CNVRT028 

CNVRT030 

CNVRT032 

CNVRT034 

CNVRT036 

GNVRT038 

CNVRT040 

CNVRT042 

CNVRT044 

CNVRT046 

CNVRT048 

CNVRT050 

CNVRT052 

CNVRT054 

CNVRT056 

CNVRT058 

CNVRT060 

CNVRT062 

CNVRT064 

CNVRT066 

CNVRT068 

CNVRT070 

CNVRT072 

CNVRT074 

CNVRT076 

CNVRT078 

CNVRT080 

CNVRT082 

CNVRT084 

CNVRT086 


ii 






ft 


ii 


v 

I 
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c 

c 

c 

c 

c 

c 

c 


*** MAIN FOR INVARA USING ALLMAG AND THE COORDINATE CONVERSION 


*** DOUBLE PRECISION DECK FOR IBM 360 MACHINES (EBCDIC, 029 PUNCH) 
*** INPUT: 


*** 

*** 

*** 

*** 


1 


2 

3 

4 

5 


7 


8 


9 

99 


ICOORD REFERENCE SYSTEM OF INPUT COORDINATES 

=1 GEODETIC; =2 GEOCENTRIC. 

MODEL CHOICE OF 7 MODELS (FROM ALLMAG) 

TM TIME IN YEARS FOR DESIRED FIELD 

IMPLICIT REAL*8 (A-H,0-Z) . 

DATA RAD/57.2957795/ 

FORMAT( 'O' ,9X, •GEOCENTRIC COORDINATES GEODETIC COORDINATES 1 , 

1 10X, 1 B * ,14X» 1 L # , / , I IX , ' L AT LONG ALTS7X,'LAT LONGS 

2 4X , 1 ALT S / , 10X , 1 ( DEGR ) ( DEGR ) ( KM ) 1 ,6X , 1 ( DEGR ) ( DEGR ) (KM)S 

3 9X, 1 (GAUSS) S6X, * (EARTH RADII)*) 

FORMAT ( 215, FLO. 3 ) 

FORMAT ( •0 , ,7X,2F8.2,F8.1,1X,2F8.2,F8.1,2X,F12.5,F12.4) 
F0RMAT(3F10.3) 

FORMAT ( 1 1 MODEL S 13,/, * T I ME * , F8. 1 , / , • ICOORD', 13) 

READ (5,2) ICOORD, MODEL, TM 
WR I TE ( 6, 5 ) MODEL, TM, ICOORD 
WR I TE ( 6 , 1 ) 

I F ( ICOORD. EQ. 2) GO TO 8 

READ( 5,4,END=99) GDL AT , GLON ,GDALT 

CALL CONVRTI 1 , GDL AT , GDALT , GCL AT , RKM ) 

CT=DSIN(GCLAT/RAD) 

GCALT=RKM-6378.16/DSQRT ( 1 . +. 0067397*CT*CT ) 

GO TO 9 

READ ( 5,4, END=99) GCL AT , GLON , GC ALT 
CT=DSIN(GCLAT/RAD) 

RKM=GCALT+6378. 16/DSQRT ( 1 • +• 0067397*CT*CT ) 

CALL CONVRTt 2, GDL AT , GDALT , GCL AT, RKM ) 

CALL INVARA(M0DEL,TM,GCLAT,GL0N,GCALT,0.01D0, BB,FL ) 

WRITE (6, 3 ) GCL AT, GLON, GC ALT, GDL AT, GLON, GDALT, BB,FL 

GO TO (7,8), ICOORD 

STOP 

END 


002 

004 

006 

008 

010 

012 

014 

016 

018 

020 

022 

024 

026 

028 

030 

032 

034 

036 

038 

040 

042 

044 

046 

048 

050 

052 

054 

056 

058 

060 

062 

064 

066 

068 

070 


2 1 

50.000 

30.000 

10.000 
- 10.000 
-30.000 
-50.000 

50.000 

30.000 

10.000 
- 10.000 
-30.000 
-50.000 

50.000 

30.000 

10.000 
- 10.000 
-30.000 
-50.000 

50.000 

30.000 
- 10.000 
-30.000 
-50.000 

10.000 


1965.0 

60.000 

60.000 

60.000 

60.000 

60.000 

60.000 

160.000 

160.000 

160.000 

160.000 

160.000 

160.000 

260.000 

260.000 

260.000 

260.000 

260.000 

260.000 

360.000 

360.000 

360.000 

360.000 

360.000 

360.000 


1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 

1000.000 


Input data used in INVARA 
test run, chosen to be 
compatible with Table 11-2 
in Hassitt & Mcllwain (1967) 


NASA-GSFC COML., Arlington, Va. 
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OUTPUT FROM STANDARD CALLING ROUTINE FOR GDALMG 


MOO CL 

TIME 

LAT 

LONG 

ALT 

X 

Y 

z 

F 

M 

DEC 

INC 

MODEL 

USED IS NUMBER I 

HENDR ICKSCCAIN 

99-TEPM GSFC 9/65 FOR TM - 1965 

i.ooo 




t 

1963, 

-30. 

-90. 

100. 

0.23810 

0.06764 

-0.16808 

0.29919 

0.24752 

13*86 

-34.10 

k 

1965. 

-30. 

90. 

100. 

0.19959 

-0.06856 

-0.45810 

0.50437 

0.21 104 

-18*96 

-65.27 

1 

1965. 

60. 

-90. 

100. 

0.05993 

-0.00630 

0.57745 

0.58059 

0.06026 

-6.01 

84*04 

t 

1965. 

60 . 

90. 

100. 

0. 12668 

0.0211S 

0.56277 

0.57724 

0.12843 

9*48 

77.14 

MOO EL 

used is NUMBER 1 

HENDR ICKSG CAIN 

'99-TERN GSFC 9/65 FOR TM « 1970 

i.OOO 




t 

1970. 

-30. 

-90. 

100. 

0.23523 

0.06664 

-0.16583 

0.29543 

0*24448 

15*82 

-34.15 

1 

1970. 

-30. 

90. 

100. 

0. 19883 

-0. 06870 

-0.46040 

0.50619 

0*21038 

-19*06 

-63*44 

1 

1970. 

60. 

-90. 

100. 

0.06146 

-0.00600 

0.57795 

.0.58124 

0*06176 

-5*57 

83.90 

! 

1970. 

60. 

90. 

100. 

0.12698 

0.02062 

0.56456 

0.57903 

0* 12864 

9.22 

77.16 

MODEL 

USED IS 1 

NUMBER 2 

CAIN ET 

•AL. 120-TERM GSFC 

12/66 FOR TM * 1965.000 




Z 

1965. 

-30. 

-90. 

100* 

0.23823 

0.06681 

-0.16840 

0.29929 

0.24742 

15.67 

-34,24 

z 

1965. 

-30. 

90. 

100. 

0. 20058 

-0.06847 

-0.45654 

0.50334 

0.21 194 

-18*65 

-65*10 

2 

1965. 

60. 

-90. 

100. 

0.06005 

-0.00627 

0.57762 

0.58077 

0.06036 

-5.96 

64.03 

2 ‘ 

1965. 

60. 

90. 

100. 

0.1271 6 

0*01965 

0.56235 

0*57688 

0. 12667 

8*78 

77.11 

MOO EL 

USED IS 1 

NUMBER 2 

CAIN ET.AL. 120-TERM GSFC 

12/66 FOR TM - 1970.000 




2 

1970. 

-30. 

-90. 

100. 

0.23521 

0.06517 

-0.16674 

0 .29559 

0.24407 

15*49 

-34*34 

2 

1970. 

-30. 

90. 

100. 

0.20194 

-0.06809 

-0.45721 

0 • 50444 

0*21311 

-18.63 

-65.01 

2 

1970. 

60. 

-90. 

100. 

0.06205 

-0. 00579 

0.57831 

0.58166 

0*06232 

-5.33 

83.85 

2 

1970. 

60. 

90. 

100. 

0* 1 2795 

0.01917 

0.56524 

0.57986 

0. 12938 

8*52 

77.11 

MODEL 

USED IS NUMBER 3 

CAINCLANGEL 143-TERM POGO 

10/68 FOR TM - 1965.000 




3 

1965. 

-30. 

-90. 

100. 

0.2 395 7 

0.06703 

-0. 16721 

0.29975 

0*24877 

15*63 

-33.91 

3 

1965. 

-30. 

90. 

100. 

0.20044 

-0.06840 

-0.45674 

0. SO 3 46 

0*21179 

— 16 * 8 A 

-65*12 

3 

1969. 

60 . 

-90. 

100. 

0.0 5991 

-0.00639 

0.57780 

0.58093 

0*06025 

-6.09 

84.05 

3 

1965. 

60. 

90. 

100. 

0. 12700 

0.01963 

0.56318 

0*57766 

0* 12850 

8*79 

77.15 

MOOEL 

USED IS 1 

NUMBER 3 

CAINCL ANGEL 1 

43- TERM POGO 

t 10/68 FOR TM * 1970.000 




3 

1970. 

-30. 

-90. 

100. 

0. 23700 

0.06798 

-0.16310 

0*29562 

0*24655 

16*00 

-33.49 

3 

1970. 

-30. 

90. 

100. 

0. 19989 

-0.06857 

-0.45668 

0*50321 

0*21132 

-18*93 

-63*17 

3 

1970. 

60. 

-90. 

100. 

0.06154 

-0.00575 

0.57938 

0.58267 

0*06101 

-5.33 

83.91 

3 

1970. 

60. 

90. 

100. 

0. 1 2860 

0*01883 

0.56270 

0.57751 

0.12997 

8*33 

76.99 

MOOEL 

USED IS NUMBER 4 

CAINCSMEENEV 

120-TERM POGO 8/69 FOR TM - 1965.000 




♦ 

1963. 

-30. 

-90. 

100. 

0.23954 

0.06687 

r 0. 16707 

0.29960 

0.24870 

15.60 

-33.89 

4 

1963. 

-30. 

90. 

100. 

0.20046 

-0.06802 

—0 .45679 

0.50346 

0.21168 

-18.74 

-65*14 

4 

1965. 

60. 

-90. 

100. 

0.05987 

-0.00635 

0.57778 

0.58091 

0.06021 

-6.05 

84*05 

4 

1965. 

60. 

90. 

100. 

0. 12692 

0.01966 

0.56312 

0*57759 

0. 12844 

8*80 

77.13 

MOOEL 

USED IS 1 

NUMBER 4 

CA IN LS WEE NET 

120-TERM POGO 8/69 FOR TM ■ 1970.000 




4 

1970. 

-30. 

-90. 

100. 

0.23596 

0.06464 

-0.16580 

0*29555 

0.24466 

13*32 

-34.13 

4 

1970. 

-30. 

90. 

100. 

0.20022 

-0.06888 

-0.45658 

0*50329 

0.21174 

-18*98 

-65*12 

4 

1970. 

60. 

-90. 

too. 

0.06144 

-0.00616 

0.57940 

0.58269 

0.06175 

-5.73 

63*92 

4 

1970. 

60. 

90. 

100. 

0. 1 2 833 

0.01870 

0.56263 

0.57738 

0.12969 

6*29 

77.02 

MODEL 

USED IS 

NUMBER 3 

1 GRF 

1965. 

0 80-TERM 10/68 FOR TM « 1965.000 




3 

1965. 

-30. 

-90. 

100. 

0.23816 

0.06658 

-0.16682 

0.29830 

0.24730 

15*62 

-34.00 

3 

1965. 

-30. 

90. 

100. 

0. 20078 

-0.06798 

-0.45695 

0.50372 

0.21 198 

-18.71 

-65*1 1 

3 

1965. 

60. 

-90. 

100. 

0.05967 

-0.00574 

0.57735 

0.58045 

0.05995 

-5*50 

84.07 

3 

1963. 

60. 

90. 

100. 

0. 1 2650 

0.01888 

0.S6304 

0.57739 

0.12790 

8*49 

77.20 

MODEL 

USED IS 

KWMBER 5 

IGRF 

1965. 

0 80-TERM 10/68 FOR TM » 1970.000 




3 

1970. 

-30. 

-90. 

100. 

0.2 3596 

0.06537 

-0.16451 

0.29498 

0.24484 

15*48 

-33.90 

3 

1970. 

-30. 

90. 

100. 

0. 20077 

-0.06783 

-0.45890 

0.50547 

0*21 192 

-18*67 

-65.21 

3 

1970. 

60. 

-90 • 

100. 

0.06126 

-0.00540 

0.57805 

0.58132 

0*06150 

-5.04 

83.93 

5 

1970. 

60. 

90. 

100. 

0*12723 

0.01825 

0.56485 

0.57929 

0* 12853 

8* 16 

77.18 

MOOEL 

USED IS 

NUMBER 6 

LEAToN 

MALIN 

EVANS 80-TERM 1965 FOR TM * 1965.000 




6 

1965. 

-30. 

-90. 

100. 

0.24028 

0.06833 

-0.16703 

0.30050 

0*24980 

15*87 

-33.77 

6 

1963. 

-30. 

90. 

100. 

0.20050 

-0.06946 

-0.45780 

0.50459 

0*21219 

-19*11 

-65.13 

6 

1963. 

60. 

-90. 

10O. 

0.05998 

-0.00619 

0*57939 

0.58251 

0*06030 

-5.89 

84.06 

6 

1963. 

60. 

90. 

100. 

0. 12554 

0.01928 

0.56330 

0.57744 

0*12701 

8.73 

77.29 

MOOEL 

USED IS 

NUMBER 6 

LEATON 

MAL IN 

EVANS 80- TERM 1963 FOR TM * 1970.000 




e 

1970. 

-30. 

-90. 

100. 

0*2 3841 

0.06697 

-0.16406 

0.29706 

0*24764 

13.69 

-33.52 

6 

1970. 

-30. 

90. 

100. 

0* 19879 

-0.06963 

-0.46038 

0 . 506 28 

0*21064 

-19*31 

—65. At 

6 

1970. 

60. 

-90. 

lOO. 

0.06155 

-0.00593 

0.50080 

0*58408 

0*06183 

-5*51 

83.92 

• 

1970. 

60. 

90. 

100. 

0. 12651 

0.01856 

0.56469 

0.57899 

0* 12786 

8*34 

77.24 

MOOEL 

USED IS 

NUMBER 7 

HURKITZ US 

CtGS 1 60- TERM 1970 FOR TM - 1965.000 




7 

1965. 

-30. 

-90. 

100. 

0.23904 

0.06707 

-0.16777 

0.29965 

0*24828 

15.67 

-34.05 

7 

1965. 

-30. 

90. 

100. 

O. 1994B 

-0.06660 

-0.45625 

0.50266 

0*21094 

-18*96 

—65. 1 9 

7 

1965. 

60. 

-90. 

100. 

0.06025 

-0.0061 7 

0.57719 

0.58036 

0.06057 

-5.85 

84.01 

7 

1965. 

60. 

90. 

100. 

0. 1 2699 

0.01995 

0.56219 

0.57670 

0*12855 

8*03 

77.12 

MOOEL 

USED IS 

NUMBER 7 

HURWITZ US 

CSGS 168- TERM 1970 FOR TM » 1970.000 




7 

1970. 

-30. 

-90. 

100. 

0.23630 

0.06532 

-0.16498 

0.29550 

0.24516 

15*45 

-33.94 

7 

1970. 

-30. 

90. 

100. 

0. 19896 

-0.06348 

-0.45681 

0.50294 

0.21042 

-18.99 

-65.27 

7 

1970. 

60. 

-90. 

100. 

0.06184 

-0.00597 

0.57796 

0.58129 

0.06213 

-5.51 

83.86 

7 

1970. 

60. 

90. 

lOO. 

0* 12816 

0.01953 

0.56247 

0.57721 

0.12964 

8*66 

77.02 
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1 965.00 


OUTPUT FROM LINTRA TEST RUN #1 (iPRlNT - 1} 


MOO EL 
TIME 
PRINT 
COORD 


I 


GEOCENTRIC COORDINATES 
LAT LONG RKM 

(OEGR) (DEGR) (KM) 


ORIGIN 44.8* >90.00 6367.5 


GEODETIC COORDINATES STEPS IZ E/ARCLENGTM DIR HALT LABB. 

LAT LONG ALT DS/ARC 

(OEGR) ( OEGR) (KM) (KM) 

45.00 -90.00 0.0 SIS. I. 0. TEST I 

STEP LAT LON RKM ALT BR BT 


BP 


B 



MODEL U5E0 IS NUMBER 1 HENORI CKStCAtN 99- TERM GSFC 9/65 FOR TM * 1965.000 

OLD COORDINATES FOR STEP# I • * 44. BOB -90.000 636S. 

OLD COORDINATES FOR STEP# 2 *• 44. BOS -90.000 636S. 

OLD COORDINATES FOR STEP# 3 •• 44. BOS -90.000 6368* 

OLD COORDINATES FOR STEP# 4 ** 42.B6S -90.153 7145. 

OLD COORDINATES FOR STEP# 5 ** 42.865 -90.153 7145. 

OLD COORDINATES FOR STEP# 6 ** 40.897 -90.318 7913. 

OLD COORDINATES FOR STEP# 7 ** 40.B97 -90 .318 7913. 

OLD COORDINATES FOR STEP# 6 ** 38.934 -90.484 8672. 

OLD COORDINATES FOR STEP# 9 ** 36.980 -90.652 9421. 

OLD COORDINATES FOR STEP# 10 ** 35.043 -90.816 10161. 

OLD COORDINATES FOR STEP# 11 •* 33.124 -90.978 10891. 

OLD COORDINATES FOR STEP# 12 ** 31.218 -91.138 116Q9. 

OLD COORDINATES FOR STEP# 13 ** 29.323 -91.297 12316. 

OLO COORDINATES FOR STEP# 14 •* 27.437 -91.455 13011. 

OLD COORDINATES FOR STEP* 15 ** 25.555 -91.613 13692. 

OLO COORDINATES FOR STEP# 16 *» 23.674 -91.770 14358. 

OLO COORDINATES FOR STEP# 17 *# 21.791 -91.929 15008. 

OLD COORDINATES FOR STEP# 18 ** 19.901 -92.088 15641. 

OLO COORDINATES FOR STEP# 19 • * 18.003 -92.249 16254. 

OLO COORDINATES FOR STEP# 20 • • 16.092 -92.411 16846. 

OLO COORDINATES FOR STEP# 21 ** 14.166 -92.576 17414. 

OLO COORDINATES FOR STEP# 22 *m 12.222 -92.743 17957. 

OLD COORDINATES FOR STEP# 23 ** 10.257 -92.914 18471. 

OLO COORDINATES FOR STEP# 24 *• 8.270 -93.088 18953. 

OLO COORDINATES FOR STEP# 25 *# 6.258 -93.267 19400. 

OLO COORDINATES FOR STEP# 26 ** 4.221 -93.450 19809. 

OLD COORDINATES FOR STEP# 27 •• 2.158 -93.639 20175. 

OLD COORDINATES FOR STEP# 28 *4 0.070 -93.834 20496. 

OLO COORDINATES FOR STEP# 29 ** -2.041 -94.034 20767. 

OLO COORDINATES FOR STEP# 30 ** -4.174 -94.242 20984. 

OLD COORDINATES FOR STEP# 31 ** -6.323 -94.457 21145. 

OLO COORDINATES FOR STEP# 32 ** -8.485 -94.679 21247. 

OLD COORDINATES FOR STEP# 33 ** -10.653 -94.909 21288. 

OLD COORDINATES FOR STEP# 34 ** -12.822 -95.147 21268. 

OLO COORDINATES FOR STEP# 35 -14.985 -95.394 21187. 

OLD COORDINATES FOR STEP# 36 ** -17.136 -95.650 21046. 

OLO COORDINATES FOR STEP# 37 #* -19.270 -95.916 20848. 

OLD COORDINATES FOR STEP# 38 ** -21.383 -96.192 20396. 

OLD COORDINATES FOR STEP# 39 *4 -23.473 -96.480 20293. 

OLO COORDINATES FOR STEP# 40 44 -25.536 -96.779 19943. 

OLD COORDINATES FOR STEP# 41 44 -27.572 -97.093 19549. 

OLD COORDINATES FOR STEP# 42 44 -29.582 -97.422 19116. 

OLO COORDINATES FOR STEP# 43 44 -31.565 -97.769 18647. 

OLD COORDINATES FOR STEP* 44 4* -33.523 -98.137 18145. 

OLD COORDINATES FOR STEP# 45 44 -35.458 -98.527 17614. 

OLD COORDINATES FOR STEP* 46 44 -37.373 -98.945 17056. 

OLO COORDINATES FOR STEP# 47 *4 -39.269 -99.394 16474. 

OLO COORDINATES FOR STEP* 48 44 -41.149 -99.880 15871. 

OLO COORDINATES FOR STEP# 49 *» -43.017 -100.409 15248. 

OLD COORDINATES FOR STEP* 50 *« -44.675 -130.990 14607. 

OLO COORDINATES FOR STEP* 51 4* -46.726 -101.633 13949. 

OLO COORDINATES FOR STEP* 52 44 -46.572 -102.352 13278. 

OLO COORDINATES FOR STEP# 53 44 -50.417 -103.162 12592. 

OLD COORDINATES FOR STEP* 54 44 -52.264 -104.086 11895. 

OLD COORDINATES FOR STEP# 55 44 -54.114 -105.152 11187. 

OLD COORDINATES FOR STEP# 56 44 -55.968 -106.399 10468. 

OLD COORDINATES FOR STEP# 57 44 -57.828 -107.878 9741. 

OLO COORDINATES FOR STEP* 58 44 -59.692 -109.663 9005. 

OLD COORDINATES FOR STEP* 59 44 -61.555 -111.853 8262. 

OLD COORDINATES FOR STEP* 60 44 -63.404 -114.587 7S12. 

□LO COORDINATES FOR STEP# 61 44 -65.216 -118.053 6756. 

INTRSCT -66.13 -120.24 6360.3 -66.27 -120.24 0.0 45819. 


-0. -0.57342 -0.15686 0.00830 0.59455 

-1. -0.38375 -0.12143 0.00732 0.40257 

-1. -0.38448 -0.12173 0.00736 0.40335 

777. -0.38448 -0.12173 0.00736 0.40*35 

776. -0.27011 -0.09618 0.00616 0.2867? 

1544. -0.27016 -0*09614 0.00615 0.26662" 

1543. -0.19673-0.07728 0.00510 0.211*3 

2302. -0.19671 -0.07729 0.00510 0.2lf4l 

3051 • -0.14729 -0.06319 0.00427 0.160*3 

3790. -0.11281 —0.05247 0.00362 0.124*7 

4519. -0.08803 -0.04417 0.00311 0.09654 

5237. -0.06977 -0.03766 0.00270 0.07932 

5943. -0.05602 -0.03247 0.00237 0.06479 

6637. -0.04549 -0.0282C 0.00210 0.053T60 

7318. -0.03728 -0.02487 0.00188 0.04*85 

7983. -0.03080 -0.02206 0.00170 0.03792 

8633. -0.02562 -0*01972 0.00154 0.032*6 

9265. -0.02142 -0.01777 0.00141 0.02786 

9878. -0.01799 -0.01612 0.00130 0.02419 

10469. -0.01516 -0.01472 0.00121 0.02116 

1 tO 38 . -0.01279 -0.01353 0.00113 0.0186$~ 

1 1580. -0.01081 -0.01252 0.00106 0.01657 

12093. -0.00912 -0.01165 0.00100 0.01483 

12575. -0.00768-0.01091 0.00095 0.0133$ 

13022. -0.00643 -0.01029 0.00091 0.01216 

13431. -0.00534 -0.00975 0.00088 0.01116 

1 3797. -0.0 0 4 38 - 0.0C931 0.00086 0.0103? 

14118. -0.00352 -0.00894 0.00084 0.00*6* 

14388. -0.00275 -0.00864 O.OOC83 0.00911 
14606. -0.00203 -0.0C841 0.00083 0.00869 

14767. -0.00136 -0.00824 0.00083 0.00*40 

14869. -0.00073 -0.00814 0.00084 0.0082T 

14910. -0.00011 -0.00809 0.00086 0.00814 

14891. 0.00051 -0.00810 0.0QC88 0.00817 

14810. 0.00114 -0.00817 0.00092 0.00830 

14670. 0.00179 -0.00830 0.00096 0.0085* 

14473. 0.00248 -0.0C849 0.00102 0.0089* 

14221. 0.00322 -0.0C874 0.00109 0.0093* 

13918. 0.00403 -0.00908 0.00117 0.00*98 

13569. 0.00493 -0.00945 0.00127 0.01074 

1 3176. 0.00595 -0.0C992 0.00140 0.01165 

12743. 0.00711 -0.01047 0.00154 0.01275 

12275. 0.00843 -0.01112 0.00172 0.01407 

11774. 0.00998-0.01189 0.00193 0.01564 

11243. 0.01178-0.01279 0.00218 0.01752 

10686. 0.01391 -0.01383 0.00250 0.0X977 

10105. 0.01644 -0.01505 0.00288 0.02248* 

9502. 0.01949 -0.0!64e 0.00336 0.02S7S 

8879. 0.02319 -0.01817 3.00395 O.03972T 

8239. 0.02771 -0.02017 0.00471 0.03*59 

7583. 0.03331 -0.02255 0.00567 0.04062 

6911. 0.04032 -0.02540 0.00694 0.04815 

6227. 0.04920 -0.02886 0.00860 0.0576$ 

5530. 0.06064 -0.03309 0.0108S 0.06993 

4823. 0.07260 -0.03831 0.01395 0.08589 

4105. 0.09555 -0.04484 0.01830 0.10712 

3378. 0.12275 -0.0S311 0.02461 O.X36O0 

2643. 0.16088 -0.06372 0.03400 0.17635 

1900. 0.21611 -0.07753 0.04849 0.23465 

1151. 0.29962 -0.09570 0.07166 0.32259 

396. 0.43309 -0.11978 0.11018 0.46266 




i 


1 

1 

MODEL 1 

TIME 1965.0 


OUTPUT FROM 

INVARA TEST RUN 



1 

ICOORD 2 








2 

GEOCENTRIC COORDINATES 

GEODETIC COORDINATES 

* 

B 

* 

L 

3 

B 

LAT 

LONG 

ALT 

LAT 

LONG 

ALT 



l 

(DEGR) 

(DEGR) 

(KM) 

(DEGR) (DEGR) 

( KM) 

(GAUSS ) 

(EARTH F 

a 

( 

MODEL USED IS NUMBER 1 

HENDRICKSeCAtN 99 

'-TERM GSFC 9/65 FOR 

TM = 1965. 

000 

at 

50.00 

60.00 

1000.0 

50.16 

o 

o 

t 

o 

1000.0 

0.34101 

2.2417 


30.00 

60.00 

1000.0 

30.14 

60.00 

1000.0 

0.28446 

1.3330 

• ' 

10.00 

60.00 

1000.0 

10.06 

60.00 

1000.0 

0.23211 

1 .1092 


-10.00 

60.00 

1000.0 

-10.06 

60.00 

1000.0 

0.23609 

1 .2581 

- 

-30.00 

60.00 

1000.0 

-30.14 

60. 00 

1000.0 

0.26279 

1.9259 

- 

-50.00 

60.00 

1000.0 

-50.16 

60.00 

1000.0 

0.29988 

3.9519 

a 

50.00 

160.00 

1000.0 

50. 16 

160.00 

1000.0 

0.32219 

2.1171 


30.00 

160.00 

1000.0 

30.14 

160.00 

1000.0 

0.25390 

1 .3242 

- 

10.00 

160.00 

1000.0 

10.06 

160. 00 

1000.0 

0.22232 

1 .1209 

.1 

K 

a 

-10.00 

160.00 

1000.0 

-10.06 

160.00 

1000.0 

0.26375 

1.2316 

1 

-30.00 

160.00. 

1000.0 

-30. 14 

160.00 

1000.0 

0.34381 

1 .8405 

3 

-50.00 

160.00 

1000.0 

-50.16 

160.00 

1000.0 

0.40628 

4.3908 


50.00 

260.00 

1000.0 

50.16 

260.00 

1000.0 

0.37821 

4.6888 


30.00 

260.00 

1000.0 

30. 14 

260.00 

1000.0 

0.32013 

1 .9539 


10.00 

260.00 

1000.0 

10.06 

260. 00 

1000.0 

0.24044 

1.2910 

- 

-10.00 

260.00 

1000.0 

-10.06 

260. 00 

1000.0 

0. 19753 

1 .1643 


-30.00 

260.00 

1000.0 

-30.14 

260. 00 

1000.0 

0.22004 

1.3276 


-50.00 

260.00 

1000.0 

-50.16 

260.00 

1000.0 

0.27878 

1 .9094 

n 

1 

50.00 

360.00 

1000.0 

50.16 

360.00 

1000.0 

0.31112 

2.6404 

1 

L 

30.00 

360.00 

1000.0 

30.14 

360. 00 

1000.0 

0.25265 

1.4280 


10.00 

360.00 

1000.0 

10.06 

360.00 

1000.0 

0.19930 

1.1712 


-10.00 

360.00 

1000.0 

-10.06 

360.00 

1000.0 

0. 18892 

1 .2951 


-30.00 

360.00 

1 000.0 

-30.14 

360.00 

1000.0 

0.19635 

1 .7096 


-50.00 

360.00 

1000.0 

-50.16 

360.00 

1000.0 

0.22230 

2.5446 


‘Compare these values of B and L with Table 11-2 in Hassitt and Mcllwain (1967) 


L-20 


